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1.0  ABSTRACT 


This  report  describes  a method  for  calculating  viscous  effects  on  the  lifts 

✓ 

and  pressure  distributions  of  arbitrary  three-dimensional  configurations.  The 
approach  consists  of  combining  a panel  method,  which  calculates  potential  flow 
about  arbitrary  three-dimensional  lifting  configurations,  with  a boundary-layer 
method.  Combined  procedures  have  been  constructed  using  a two-dimensional 
boundary-layer  method  in  a strip-theory  sense  and  using  a three-dimensional 
small  cross-flow  boundary-layer  method.  Various  fundamental  and  procedural 
aspects  of  the  general  calculation  scheme  are  investigated  and  discussed.  Final 
emphasis  is  on  a completely  automatic  procedure  that  performs  its  calculations 
in  a single  computer  run  without  intermediate  human  intervention.  For  this 
purpose  the  method  based  on  a strip-theory  boundary  layer  has  proved  very  satis- 
factory. Calculated  inviscid  and  viscous  lift  and  pressure  distributions  are 
compared  with  experimental  data  for  a variety  of  wings  and  wing-fuselages  having 
both  conventional  and  supercritical  airfoil  sections.  The  agreement  of  the  cal- 
culations with  experiment  appears  to  be  quite  good. 


1 


'•l  ' — 1 - v-1*'  «• 


2.0  TABLE  OF  CONTENTS 


Page 

1 .0  Abstract 1 

2.0  Table  of  Contents 2 

3.0  Index  of  Figures 5 

4.0  Principal  Notation  8 

5.0  Introduction  9 

6.0  General  Aspects  of  Procedures  for  Approximating  Viscous  Flows 

by  Means  of  Combined  Inviscid-Flow  Boundary-Layer  Methods  ....  12 

6.1  Two-Dimensional  and  Three-Dimensional  Flows.  Boundary- 

Layer  Methods 12 

6.2  Methods  of  Boundary-Layer  Simulation  13 

6.3  The  Wake 14 

6.4  Sequential  Use  of  the  Potential -Flow  and  Boundary-Layer 

Methods 14 

6.5  Comparison  of  Computational  Aspects  of  the  Two  Boundary- 

Layer  Simulation  Methods  16 

6.6  Application  of  the  Kutta  Condition 18 

6.7  Form  of  the  Surface  Vorticity  19 

7.0  Two-Dimensional  Studies  21 

7.1  General  Remarks 21 

7.2  Point-Number  Requirements  for  the  Boundary-Layer  Program  . . 21 

7.3  Effect  of  the  Wake 22 

7.4  Results  for  Selected  Airfoils  23 

7.4.1  RAE  101 23 

7.4.2  64A010 24 

7.4.3  Symmetric  Joukowsky  24 

7.5  Conclusions 25 


2 


Page 


f 


i 

8.0  A Three-Dimensional  Procedure  Based  on  a Strip-Theory  Boundary- 

Layer  Method  n 

8.1  General  Remarks 27 

8.2  Some  Aspects  of  the  Three-Dimensional  Potential -Flow  Panel 

Method 28 

8.3  General  Description  of  the  Combined  Potential -Flow 

Boundary-Layer  Method  29 

8.3.1  The  Initial  Potential -Flow  and  Boundary-Layer 

Calculation 29 

8.3.2  Surface  Displacement  Simulation  32 

8.3.3  Surface  Blowing  Simulation  33 

8.4  Comparison  of  Calculated  Results  for  Isolated  Wings  with 

Experimental  Data 34 

8.4.1  Symmetric  Swept  Wing 35 

8.4.2  Cambered-Twisted  Swept  Wing  36 

8.5  Comparison  of  Calculated  Results  for  Wing-Fuselages  with 

Experimental  Data 36 

8.5.1  A Straight  Wing  on  a Round  Fuselage 37 

8.5.2  Conventional  Swept  Wing  Mounted  Low  on  a Fuselage  ...  38 

8.5.3  Supercritical  Swept  Wing  Mounted  High  on  a 

Fuselage 39 

8.5.4  A Supercritical  Transport  Configuration  40 

8.6  Conclusions 41 

9.0  A Three-Dimensional  Procedure  Based  on  a Smal 1 -Cross-Flow 

Boundary-Layer  Method  43 

9.1  General  Remarks 43 

9.2  Surface  Coordinates  43 

I 9.3  Surface  Velocity  Components  46 

9.4  Streamline  Computation  47 

9.5  Curvature  Calculation  48 


3 


L 


r 


Page 

48 


9.6  Boundary-Layer  Calculation  and  Simulation 

9.7  Calculated  Results  

10.0  References  


IV 


3.0  INDEX  OF  FIGURES 


No.  Page 

1 Representation  of  a three-dimensional  lifting  configuration . ...  53 

2 Two  methods  of  boundary- layer  simulation  54 

3 Application  of  the  Kutta  condition  54 

4 Calculated  upper-surface  displacement-thickness  distributions  on 

an  airfoil  obtained  using  various  point  numbers  55 

5 Effect  of  the  wake  displacement  thickness  on  calculated  lift 
coefficients  for  a 11.8%-thick  Joukowsky  airfoil  at  6°  angle 

of  attack 56 

6 Comparison  of  calculated  inviscid  pressure  distributions  and  lifts 
with  exact  values  for  an  RAE  101  airfoil  at  8.2°  angle  of  attack. 

52  surface  panels 57 

7 Comparison  of  calculated  and  experimental  pressure  distributions 

and  lifts  for  an  RAE  101  airfoil  at  8.2°  angle  of  attack.  Constant 
surface  vorticity  distribution  58 

8 Comparison  of  calculated  and  experimental  pressure  distributions 

and  lifts  for  an  RAE  101  airfoil  at  8.2°  angle  of  attack.  Para- 
bolic surface  vorticity  distribution  59 

9 Comparison  of  calculated  inviscid  pressure  distributions  and  lifts 
with  exact  values  for  a 64A010  airfoil  at  8°  angle  of  attack. 

40  surface  panels  60 

10  Comparison  of  calculated  and  experimental  pressure  distributions 

and  lifts  for  a 64A010  airfoil  at  8°  angle  of  attack.  Constant 
surface  vorticity  distribution  61 

11  Comparison  of  calculated  inviscid  pressure  distributions  and  lifts 
with  exact  values  for  a Joukowsky  airfoil  at  6°  angle  of  attack. 

54  surface  panels 62 

12  Comparison  of  calculated  and  experimental  pressure  distributions 

and  lifts  for  a Joukowsky  airfoil  at  6°  angle  of  attack.  Constant 
surface  vorticity  distribution  63 

13  Comparison  of  calculated  and  experimental  pressure  distributions 

and  lifts  for  a Joukowsky  airfoil  at  6°  angle  of  attack.  Parabolic 
surface  vorticity  distribution  64 

14  Addition  of  displacement  thickness  to  a body 65 


15  Two  isolated  swept  wings 


5 


• vAP-— »H*:r  • 


rr-  » ~ 

No. 

1 

Page 

16 

Comparison  of  calculated  and  experimental  spanwise  distributions 
of  section  lift  coefficient  for  a swept  wing  with  symmetric 
airfoil  section  at  8.2°  angle  of  attack  

. 67 

17 

Comparison  of  calculated  and  experimental  chordwise  pressure 
distributions  at  55.5%  semi-span  for  a sweat  winq  with  symmetric 
airfoil  section  at  8.2°  angle  of  attack  

. 68 

18 

Comparison  of  calculated  and  experimental  spanwise  distributions 
of  section  lift  coefficient  for  a cambered  twisted  swept  wing  at 
8.2°  angle  of  attack  

. 69  1 

19 

Comparison  of  calculated  and  experimental  chordwise  pressure 
distributions  at  55.5%  semi -span  for  a cambered  twisted  swept 
wing  at  8.2°  angle  of  attack  

. 70 

20 

A straight  wing  on  a round  fuselage  

• 71 

21 

Comparison  of  calculated  and  experimental  spanwise  distributions 
of  section  lift  coefficient  for  a straight  wing  on  a round 
fuselage  at  6°  angle  of  attack  

■ 72  : 

22 

Comparison  of  calculated  and  experimental  chordwise  pressure 
distributions  at  60%  semi -span  on  a straight  wing  mounted  on  a 
round  fuselage  at  6°  angle  of  attack  

• 73 

23 

A conventional  swept  wing  mounted  as  a low  wing  on  a fuselage  . . 

• 74  i 

24 

Comparison  of  calculated  and  experimental  spanwise  distributions 
of  section  lift  coefficient  for  a conventional  swept  wing  mounted 
low  on  a fuselage  at  6.9°  angle  of  attack  

i 

. 75 

25 

Comparison  of  calculated  and  experimental  chordwise  oressure 
distributions  at  60%  semi -span  on  a conventional  swept  wing 
mounted  low  on  a fuselage  at  6.9°  angle  of  attack  

. 76 

26 

A supercritical  swept  wing  mounted  as  a high  wing  on  a fuselage  . 

. 77 

27 

Comparison  of  calculated  inviscid  pressure  distributions  and  lifts 
with  exact  values  for  a two-dimensional  supercritical  airfoil  at 
7.02°  angle  of  attack  

. 78 

28 

Comparison  of  calculated  and  experimental  spanwise  distributions 
of  section  lift  coefficient  for  a supercritical  swept  wing  mounted 
high  on  a fuselage  at  7.02°  angle  of  attack  

. 79 

29 

! 

Comparison  of  calculated  and  experimental  chordwise  pressure 
distributions  at  60%  semi -span  on  a supercritical  swept  wing 
mounted  low  on  a fuselage  at  7.02°  angle  of  attack  

. 80 

30 

A supercritical  transport  configuration  

. 81  1 

No. 


Page 


31  Comparison  of  calculated  inviscid  pressure  distributions  and  lifts 

with  exact  values  for  a three-dimensional  transport-type  super- 
critical airfoil  at  2.06°  angle  of  attack  82 

32  Comparison  of  calculated  and  experimental  spanwise  distributions  of 

section  lift  coefficient  on  the  wing  for  a supercritical  transport 
configuration  at  2.06°  angle  of  attack  83 

33  Comparison  of  calculated  and  experimental  chordwise  pressure 

distributions  at  the  50%  semi-span  location  on  the  wing  for  a super- 
critical transport  configuration  at  2.06°  angle  of  attack  84 

34  Definition  of  surface  coordinates  for  a wing 85 

35  A wing  in  its  coordinate  space 85 

36  An  individual  surface  panel  86 

37  Block  diagram  of  the  combined  potential -flow  smal  1 -cross -flow- 

boundary-iayer  program  86 

38  Calculated  streamlines  on  a swept  wing  with  symmetrical  airfoil 

section  at  8.2°  angle  of  attack 87 

39  Calculated  spanwise  distributions  of  displacement  thickness 

along  the  upper-surface  trailing  edge  of  a swept  wing  with 
synmetrical  airfoil  section  at  8.2°  angle  of  attack  88 

40  Comparison  of  calculated  and  experimental  spanwise  distributions  of 

section  lift  coefficient  for  a swept  wing  with  symmetric  airfoil 
section  at  8.2°  angle  of  attack 89 


! 


7 


4.0  PRINCIPAL  NOTATION 


C.  total  lift  coefficient:  equals  total  integrated  pressure  force  on 

a body  in  a direction  perpendicular  to  the  freestream  divided  by 
the  product  of  a reference  area  and  freestream  dynamic  pressure 

Cp  pressure  coefficient:  equals  difference  of  local  static  pressure 

from  freestream  static  pressure  divided  by  freestream  dynamic 
pressure 

c section  lift  coefficient:  equals  the  component  perpendicular  to  the 

freestream  of  the  integrated  pressure  force  on  a strip  of  elements 
on  a wing  divided  by  the  product  of  freestream  dynamic  pressure  and 
the  projection  of  the  area  of  the  strip  into  the  plane  containing 
the  chord  line  of  the  wing 

L total  arc-length  of  an  N-line 

N-line  a curve  on  the  surface  of  a body  defined  by  input  points  (figure  1) 
ri  unit  normal  vector  to  the  body  surface 

Re  Reynolds  number 

t integration  variable  for  calculating  surface  streamlines 

U freestream  velocity 

u,v  nonorthogonal  coordinates  in  a body  surface 

u^,v^  values  of  u and  v pertaining  to  the  midpoint  of  a panel 

V surface  velocity 

w width  of  a lifting  strip  of  elements  (figures  34  and  36) 

x/c  distance  along  an  airfoil  as  a fraction  of  its  chord 

x,y,z  Cartesian  coordinates 

5*  boundary-layer  displacement  thickness 

C,n  orthogonal  coordinates  in  the  plane  of  a panel  (figure  36) 
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5.0  INTRODUCTION 


At  the  present  time  powerful  and  accurate  computer  methods  exist  for 
calculating  inviscid  flow  about  two-  and  three-dimensional  configurations. 

One  class  of  methods  is  aoproximate  in  the  sense  that  some  smal 1 -perturbation 
assumptions  have  been  incorporated,  and  such  methods  are  apolicable  only  to 
a limited  class  of  geometric  configurations.  A second  class  of  methods  is 
"exact"  in  the  sense  that  no  restrictive  assumptions  are  employed,  and  such 
methods  are  applicable  to  all  configurations.  Especially  powerful  are  the 
so-called  "panel  methods"  for  low-speed  flow  [1].  This  last  class  of  methods 
takes  advantage  of  the  fact  that  only  in  the  low-speed  regime  is  the  govern- 
ing flow  eguation  linear,  and  it  employs  the  principle  of  superposition  of 
elementary  solutions.  The  result  is  that  panel  methods  need  consider  only 
the  surface  of  the  body  about  which  flow  is  to  be  computed,  as  opposed  to 
exact  methods  in  other  speed  regimes,  which  must  consider  the  entire  flow 
field  even  if  flow  on  the  body  is  the  only  interest.  Thus,  low-soeed  panel 
methods  are  very  efficient.  They  can  calculate  flow  about  a given  configura- 
tion in  far  less  time  than  other  exact  methods,  or,  alternatively,  for  reason- 
able computing  times  a panel  method  can  calculate  flow  about  a much  more 
complicated  configuration  than  other  exact  methods  can. 

Despite  the  neglect  of  viscosity,  surface  oressure  distributions  calcu- 
lated by  low-speed  panel  methods  agree  with  experiment  remarkably  well  in  a 
wide  variety  of  cases  [1].  There  are  only  two  important  exceptions:  (1) 

k lifting  flows,  and  (2)  flows  with  large  regions  of  separated  flow.  An  example 

of  the  first  class  is  the  flow  about  a lifting  wing  at  a moderate  angle  of 
r attack,  for  which  the  flow  is  essentially  unseparated.  Examples  of  the  second 

^ class  are  flow  about  a sphere  and  flow  about  a lifting  wing  near  maximum  lift. 

: Surface  pressures  on  a nonlifting  body  with  essentially  unseparated  flow  agree 

; very  closely  with  inviscid  surface  pressures.  The  calculation  of  fully  sepa- 

t rated  flows  is  currently  an  unsolved  problem,  despite  much  effort  by  various 

i 

I investigators.  The  subject  of  this  report  is  the  calculation  of  "essential ly- 

I unseparated"  flow  about  a body  for  which  the  boundary  layer  either  is  unsepa- 

i rated  or  else  separates  so  near  the  downstream  end  of  the  body  that  extrapola- 


tion of  the  displacement  thickness  to  the  end  of  the  body  is  justified.  In 
view  of  the  above  discussion,  only  a lifting  body  experiences  significant 
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viscous  effects,  and  the  method  to  be  presented  concentrates  on  this  case. 

The  requirement  also  has  been  imposed  that  the  method  be  fully  automatic. 

The  user  simply  specifies  the  geometry  of  the  body  and  the  flow  conditions, 
and  the  method  calculates  all  the  required  viscous  and  inviscid  flow  proper- 
ties without  further  action  on  the  part  of  the  user. 

Because  of  considerations  of  computing  time,  an  "exact"  numerical  solution 
of  the  Navier-Stokes  equations  is  not  practical  at  this  time.  Such  a solution 
would  require  far  more  computing  time  than  any  inviscid  method,  even  one  that 
considers  the  entire  flow  field,  let  alone  a relatively  fast  panel  method. 
Accordingly,  many  investigators  have  concentrated  on  methods  that  take  advan- 
tage of  the  fact  that  for  essentially  unseparated  flow  about  a body  at  a 
practical  Reynolds  number,  the  effects  of  viscosity  are  appreciable  only  in  a 
very  thin  boundary  layer  adjacent  to  the  body's  surface  and  in  a thin  wake 
downstream  of  the  body.  Such  methods  attempt  to  "patch  together"  an  inviscid 
solution  and  a boundary-layer  solution  to  represent  a '"eal  flow.  This  is 
the  approach  adopted  in  the  present  investigation.  Previous  papers  [2],  [3] 
have  presented  methods  of  this  type  applicable  to  two-dimensional  lifting 
flow.  The  present  work  is  an  attempt  to  construct  a similar  method  for  three- 
dimensional  lifting  flow.  Various  procedures  are  possible,  each  having  a 
different  degree  of  approximation,  and  several  have  been  tried. 


Many  of  the  ideas  to  be  applied  to  the  three-dimensional  case  have  direct 
two-dimensional  analogies.  Accordingly,  series  of  two-dimensional  studies 
have  been  conducted  to  deter.nine  the  effectiveness  of  various  orocedures  without 
incurring  the  high  three-dimensional  computing  costs.  The  general  conclusion 
to  be  drawn  is  that  in  view  of  all  the  approximations  being  applied  to  the 
calculation,  in  particular  those  associated  with  the  inviscid-flow  boundary- 
layer  "patching,"  the  simpler  procedures  for  representing  viscous  three- 
dimensional  flows  are  likely  to  be  as  accurate  as  more  elaborate  methods. 


Two  combined  potential -flow  boundary- layer  programs  have  been  constructed. 
The  first  combines  a three-dimensional  low-speed  panel  method  with  a two- 
dimensional  boundary- layer  method  applied  in  a "strip-theory"  sense.  The 


10 


method  based  on  a small  cross-flow  assumption.  The  various  sections  of  this 
report  describe  these  methods,  the  various  assumptions  employed,  and  the 
evidence  supporting  the  validity  of  the  assumptions.  Also  presented  are  cal- 
culated results  for  a variety  of  configurations  and  comparison  of  these  results 
with  experimental  data.  The  most  fully  developed  method  is  that  based  on  the 
two-dimensional  boundary-layer  calculation.  This  combined  program  is  fully 
automatic  and  essentially  without  "bugs."  Moreover,  its  predictions  agree 
with  experiment  very  well.  It  may  be  considered  an  operational  design  tool. 

The  more  elaborate  program  based  on  the  small  cross-flow  boundary- layer  tech- 
nique is  somewhat  more  difficult  to  run  and  may  be  said  to  require  an  expert 
user.  Moreover,  for  most  geometries,  no  substantial  gain  in  accuracy  over 
that  of  the  simpler  program  can  be  achieved.  Nevertheless,  the  second  program 
has  been  fully  documented  and  it  does  produce  answers,  including  both  boundary- 
layer  information  and  new  potential -flow  information  such  as  surface  streamlines. 
Presumably  the  second  program  will  prove  useful  in  cases  where  the  three- 
dimensional  nature  of  the  boundary  layer  is  essential. 
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6.0  GENERAL  ASPECTS  OF  PROCEDURES  FOR  APPROXIMATING  VISCOUS  FLOWS 
BY  MEANS  OF  COMBINED  INV I SC  ID-FLOW  BOUNDARY-LAYER  METHODS 


6 . 1 Two-Dimensional  and  Three-Dimensional  Flows.  Boundary-Layer  Methods. 

The  basic  procedure  for  combininq  an  inviscid-flow  method  and  a boundary- 
layer  method  is  the  same  for  two-dimensional  and  three-dimensional  flows  in  all 
speed  regimes.  The  dimensionality  and  the  speed  regime  affect  only  what  par- 
ticular inviscid-flow  method  and  what  particular  boundary-layer  method  are 
used  in  the  combined  procedure.  There  are,  of  course,  several  choices  of  both 
types  of  method  available  even  after  the  dimensionality  and  the  speed  regime 
a>"p  fixed. 

In  the  present  work  attention  is  restricted  to  the  low-speed  regime, 
where  either  the  flow  is  incompressible  or  compressibility  effects  are  small 
enough  to  be  accounted  for  by  simple  corrections  to  an  incompressible  flow 
method.  The  inviscid-flow  method  considered  is  the  well-known  "exact"  surface- 
panel  method  for  incompressible  potential  flow  [1].  Several  versions  of  this 
method  exist,  but  all  represent  the  surface  of  the  body  about  which  flow  is  to 
be  computed  by  some  combination  of  source  and  vortex  singularities  distributed 
on  surface  "panels"  that  represent  the  body  surface.  A three-dimensional 
lifting  body  represented  by  surface  panels  is  shown  in  Figure  1.  Two-  and 
three-dimensional  versions  of  this  method  exist.  They  are  direct  analogies  of 
each  other,  so  that  very  fast  two-dimensional  computer  calculations  may  be 
used  to  discover  important  properties  of  the  very  time-consuming  three- 
dimensional  procedure.  The  version  employed  in  the  present  work  is  the  source- 
panel  method  of  [1],  [4].  In  the  interest  of  brevity,  a detailed  description 
of  the  potential-flow  panel  method  is  not  given  here.  For  details,  reference 
should  be  made  to  [1],  [4].  Certain  properties  of  the  method  that  are  import- 
ant for  the  present  application  are  mentioned  as  needed  without  discussion. 

With  regard  to  boundary-layer  methods,  a two-dimensional  combined  program  [2] 
is  based  on  a finite-difference  boundary-layer  procedure  {see  Section  7.1). 

In  the  present  work  the  same  boundary-layer  method  has  been  incorporated  into  a 
combined  program  using  a strip-theory  approach.  In  the  combined  program  the 
potential-flow  calculation  dominates  the  computing  time,  especially  for  large 
panel  numbers.  (Potential  flow  computing  time  varies  as  the  square  or  the 
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cube  of  the  panel  number,  whereas  boundary-layer  computing  times  vary  as  the 
first  power.) 


Of  course  more  elaborate  three-dimensional  boundary-layer  methods  may  be 
combined  with  the  panel  method  of  [4].  A first  step  in  this  direction  was 
taken  in  the  present  work  by  using  a boundary -layer  method  based  on  an  assump- 
tion of  small  cross-flow  normal  to  the  potential -flow  streamlines.  (See 
Section  9.1 . ) 

6.2  Methods  of  Boundary-Layer  Simulation 

Once  the  inviscid  potential -flow  method  and  the  boundary-layer  method  have 
been  decided  upon,  there  still  remains  the  question  of  how  they  are  to  be  used 
to  arrive  at  an  approximation  to  high  Reynolds  number  viscous  flow,  which  is 
assumed  to  be  "essentially"  unseparated.  What  has  to  be  done  is  to  use  the 
quantities  calculated  by  the  boundary-layer  method  to  define  a modified  potential 
flow  that  is  the  desired  approximation.  This  process  may  be  thought  of  as  one 
of  simulation  of  the  boundary  layer  in  the  potential-flow  calculation. 

Four  methods  for  simulating  the  presence  of  a boundary  layer  in  both  two- 
dimensional  and  three-dimensional  potential  flow  are  contained  in  a fundamental 
paper  by  Lighthill  [5].  Two  of  these  are  much  more  corrmonly  used  than  the 
other  two  and  appear  to  have  a clearer  physical  significance.  These  are  the 
two  procedures  that  have  been  included  in  the  present  study.  They  are  illus- 
trated in  Figure  2.  In  the  first  method,  which  will  be  denoted  that  of 
surface  displacement  (Lighthill's  designation  is  flow  reduction.),  the  boundary- 
layer  displacement  thickness  is  added  to  the  original  body  in  the  direction  along 
the  local  normal  to  the  body  surface.  Thus,  a thicker  body  is  produced.  The 
potential  flow  about  this  modified  body  is  the  desired  modified  potential  flow 
that  approximates  viscous  flow.  In  the  second  method  a surface  blowing  dis- 
tribution is  defined  in  such  a way  that  the  dividing  streamline  of  the  flow  is 


the  same  as  the  modified  body  used  in  the  first  method  above.  (Lighthill's 
designation  is  equivalent  sources.)  The  potential  flow  about  the  original 
body  subject  to  this  nonzero  normal -veloci ty  boundary  condition  on  the  body 
surface  is  the  desired  modified  potential  flow  that  approximates  the  viscous 
flow.  In  both  simulation  methods  the  displacement  thickness  is  the  only  boundary- 
layer  quantity  that  affects  the  modified  potential  flow,  although  this  may 
require  some  involved  calculations  in  the  truly  three-dimensional  case. 
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6.3  The  Wake 


The  procedures  put  forward  by  Lighthill  [5]  assume  that  the  boundary-layer 
is  simulated  not  only  on  the  body  surface  but  also  in  the  wake,  i.e.,  every- 
where viscous  effects  are  important.  Thus,  in  the  surface  blowing  method  of 
section  6.2,  a suction  velocity  distribution  is  utilized  along  the  semi -infinite 
wake.  In  the  surface  displacement  method  the  thickened  modified  body  should  be 
augmented  by  a thin  semi-infinite  solid  extension  representing  the  wake  dis- 
placement thickness.  This  last  procedure  can  lead  to  numerical  instabilities 
in  the  source  panel  method  and  certainly  yields  increased  computing  time.  It 
probably  would  be  better  to  use  the  suction  simulation  of  the  wake  even  with 
the  surface  displacement  method.  In  any  case  accounting  for  viscous  effects 
in  the  wake  would  greatly  complicate  calculation  of  the  modified  potential 
flow. 

The  possible  difficulties  expressed  in  the  previous  paragraph  are  somewhat 
academic,  because  there  is  a more  serious  problem  connected  with  simulating  the 
wake.  At  present  there  are  no  accurate  procedures  available  for  calculating 
viscous  wakes  behind  lifting  airfoils,  so  there  is  no  accurate  means  for  obtain- 
ing the  wake  displacenent  thickness.  Of  course,  an  approximate  wake  displacement- 
thickness  distribution  could  be  used  for  all  bodies,  on  the  theory  that  small 
differences  in  the  wake  cannot  greatly  affect  the  flow  on  the  body  surface.  It 
is  customary,  however,  to  make  an  even  stronger  assumption,  namely  that  the  wake 
can  be  neglected  completely  without  significantly  affecting  the  lift  and  pres- 
sure distribution  on  the  body.  Published  two-dimensional  [2],  [3]  methods 
neglect  the  wake.  As  part  of  the  present  study,  calculations  were  performed 
for  a two-dimensional  airfoil  with  and  without  a wake,  on  which  the  displace- 
ment thickness  had  been  determined  experimentally  (section  7.3).  Differences 
in  the  computed  lifts  and  pressure  distributions  were  quite  small.  Accordingly, 
it  is  concluded  that  the  wake  may  be  safely  neglected  in  calculating  lift  and 
pressure  distribution  and  the  present  method  takes  advantage  of  this 
simplification. 

6.4  Sequential  Use  of  the  Potential -Flow  and  Boundary-Layer  Methods 

The  simulation  schemes  of  section  6.2  define  the  modified  potential  flow 
that  approximates  viscous  flow  in  terms  of  the  boundary-layer  displacement 
thickness.  This  last  must  be  obtained  from  a boundary-layer  method,  which 

14 


requires  as  input  the  velocity  distribution  at  the  edge  of  the  boundary  layer. 

This  velocity  distribution  is  not  the  velocity  distribution  obtained  by  an 
inviscid-flow  technique  for  flow  about  the  body  without  the  boundary  layer. 

Rather  it  is  the  velocity  distribution  corresponding  to  flow  about  the  body 
with  the  boundary  layer  already  simulated.  Thus,  the  key  flow  quantities  are 
mutually  dependent;  the  displacement  thickness  is  calculated  from  the  edge- 
velocity  distribution,  which  in  turn  is  calculated  from  the  modified  potential 
flow  which  depends  on  the  displacement  thickness. 

The  standard  remedy  for  this  dilemma  is  iteration.  To  fix  ideas  consider 
the  surface  displacement  method  of  boundary-layer  simulation.  (The  procedure 
for  the  surface-blowing  simulation  is  similar.)  This  method  can  be  used  iter- 
atively as  follows.  First  calculate  potential  flow  about  the  actual  body  shape. 
Second,  use  the  resulting  surface  velocity  distribution  as  input  to  a boundary- 
layer  method.  Third,  add  the  resulting  displacement  thickness  to  the  body  to 
obtain  a modified  body.  Fourth,  calculate  potential  flow  about  the  modified 
body,  and  iterate  the  entire  procedure  until  convergence  is  obtained.  Experi- 
ence in  two  dimensions  [2]  has  shown  that  usually  several  iterations  are 
required  for  convergence.  This  can  be  time-consuminq  (see  section  6.5)  and 
adds  an  element  of  uncertainty  to  the  calculation. 

It  is  very  tempting  to  terminate  the  above  process  after  one  iteration; 
i.e.,  obtain  displacement  thickness  from  the  inviscid  velocity  distribution 
on  the  actual  body,  add  this  to  the  body  to  obtain  t^e  modified  body,  and  use 
the  inviscid  velocity  distribution  on  this  first  modified  body  as  the  approxi- 
mation to  viscous  flow.  How  legitimate  is  this  simplification?  Or,  equivalently, 
how  closely  does  the  displacement  thickness  calculated  from  the  inviscid 
velocity  distribution  on  the  body  approximate  the  actual  displacement  thick- 
ness? Since  surface  velocities  are  not  radically  changed  by  adding  the  dis- 
placement thickness  to  the  body,  it  would  appear  that  the  above  is  a reason- 
able procedure.  The  main  objection  to  the  above  is  that  there  is  one  location 
at  which  addition  of  the  displacement  thickness  does  significantly  change  the 
surface  velocity.  If  the  aft  end  or  "trailing  edge"  of  the  body  or  airfoil  is 
not  cusped,  i.e.,  is  either  blunt  or  subtends  a nonzero  angle,  the  inviscid 
velocity  there  is  theoretically  zero.  After  the  displacement  thickness  has 
been  added  to  the  body,  the  velocity  at  the  trailing  edge  is  near  that  of  the 
free  stream  - a large  change  but  a very  local  one.  The  preceedinq  objection 


is  not  valid  for  a body  with  a cusped  trailing  edge  where  the  inviscid  veloc- 
ity is  near  free-stream  or  for  any  body  on  which  the  velocity  distribution  near 
the  trailing  edge  is  "faired  out"  to  eliminate  the  very  local  stagnation  region 
The  situation  is  unclear  for  this  case.  Van  Dyke  [6]  treats  this  problem  from 
a fundamental  standpoint  by  means  of  matched  asymptotic  expansions,  but  unfor- 
tunately only  for  the  case  of  a flat  plate  atzero  angle  of  attack.  For  this 
simple  case,  the  conclusion  is  that  iterating  the  above  boundary-layer  simula- 
tion procedure  is  not  justified. 


When  the  panel  method  of  [1],  [4]  is  used  to  calculate  flow  about  an  air- 
foil or  other  body,  velocity  is  not  calculated  at  the  trailing  edge  itself. 

The  nearest  point  to  the  trailing  edge  where  velocity  is  calculated  is  the  middle 
of  the  panel  whose  downstream  edge  is  the  trailing  edge.  In  three-dimensional  . 

cases,  computing  time  restraints  dictate  that  no  more  than  about  60  panels  be 
used  on  a wing  section,  so  that  the  midpoints  of  the  panels  adjacent  to  the 
trailing  edge  are  at  least  one-percent  of  the  airfoil  chord  away  from  the 
trailing  edge.  Because  of  the  local  nature  of  the  trai ling-edge  stagnation 
region,  the  velocity  at  these  midpoints  is  about  90  percent  of  free-stream 
velocity.  Thus  because  of  the  crudeness  of  the  numerical  representation,  all 
airfoils,  even  those  with  finite  trail ing-edge  angles,  have  velocity  distribu- 
tions similar  to  those  on  cusped  airfoils.  Therefore,  as  a practical  matter, 
iteration  of  the  boundary- layer  simulation  procedure  is  of  questionable  value 
in  calculatino  lifts  and  nressure  distributions. 


A limited  number  of  two-dimensional  cases  were  run  using  the  above  iter- 
ative procedure.  No  consistent  gain  in  accuracy  was  achieved  by  iteration 
compared  to  the  results  obtained  after  one  iteration.  Thus,  none  of  the  cases 
presented  below  continued  the  iterations  past  the  first. 


6.5  Comparison  of  Computational  Aspects  of  the  Two  Boundary-Layer  Simulation 
Methods 


From  the  discussion  of  section  6.4,  it  is  evident  that  the  procedure  for 
approximating  a viscous  flow  by  a potential  flow  with  boundary- layer  simulation 
requires  at  least  one  boundary-layer  calculation  and  two  potential -flow  calcu- 
lations: flow  about  the  actual  body  and  the  modified  potential  flow.  (If  the 

process  were  iterated  n times,  it  would  require  n boundary-layer  and 
(n  + 1)  potential -flow  calculations.)  This  is  true  for  either  method  of 


boundary-layer  simulation.  However,  the  required  computational  effort  is  not 
the  same  for  the  two  methods.  A procedure  based  on  surface  blowing  is  faster, 
because  the  second  potential -flow  calculation  (and  all  subsequent  potential- 
flow  calculations  if  iteration  is  used)  can  use  many  of  the  quantities  computed 
during  the  first  potential -flow  calculation. 

The  basic  panel-method  calculation  [1],  [4]  consists  of  two  major  parts, 
which  together  comprise  over  90  percent  of  the  required  computing  time.  The 
first  is  the  calculation  of  the  "matrix  of  influences"  that  express  the  vel- 
ocities induced  by  the  panels  on  each  other.  The  second  is  solution  of  a set 
of  linear  equations  for  the  singularity  strengths  on  the  panels.  Both  the 
"matrix  of  influences"  and  the  coefficient  matrix  of  the  linear  equations  have 
order  equal  to  the  total  number  of  panels,  and  both  matrices  are  "full,"  i.e., 
neither  contains  a large  number  of  zeros. 

The  basic  difference  between  the  two  methods  of  boundary- layer  simulation 
is  that  the  surface  blowing  method  leaves  the  geometry  of  the  body  unchanged 
and  changes  the  boundary  condition  (from  zero  to  nonzero  normal  velocity), 
while  the  surface  displacement  method  leaves  the  boundary  conditions  unchanged 
but  changes  the  geometry.  Both  the  "influence-matrix"  and  the  coefficient 
matrix  of  the  linear  equations  depend  only  on  the  geometry.  Thus,  when  calcu- 
lating the  modified  potential  flow,  the  surface  displacement  procedure  must 
recalculate  these  matrices;  while  the  surface  blowing  procedure  need  not. 
Moreover,  if  a direct  matrix  solution,  e.g.,  Gaussian  elimination,  is  used  to 
solve  the  linear  equations,  information  equivalent  to  the  inverse  of  the 
coefficient  matrix  can  be  saved  from  the  first  potential -flow  solution  and 
used  in  the  modified  potential  flow  for  the  surface  blowing  procedure  where 
the  matrices  remain  unchanged.  The  result  is  that  the  procedure  based  on 
surface  blowing  requires  only  a little  more  computing  time  than  a single  poten- 
tial flow  solution  (even  if  iteration  is  used).  On  the  contrary  the  procedure 
based  on  surface  displacement  must  perform  a completely  new  calculation  for 
the  modified  potential  flow,  and  thus  the  computing  time  is  doubled  compared 
to  a single  inviscid  case.  (If  n iterations  are  used  the  factor  is  n + 1.) 
The  above  assumes  all  computing  times,  including  that  of  the  boundary- layer 
method,  are  small  compared  to  the  potential -flow  computing  time.  This  appears 
to  be  a good  approximation.  If  an  iterative  matrix  solution  is  used  for  the 
linear  equations,  both  boundary-layer  simulation  procedures  must  solve  the 
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linear  equations  for  the  modified  potential  flow  from  the  beginning.  In  that 
case  the  computational  advantage  of  the  surface  blowing  procedure  comes  only 
from  not  having  to  form  new  matrices.  It  is  evident  that  the  surface  blowing 
method  of  boundary-layer  simulation  is  much  preferable  to  the  surface  dis- 
placement method  from  the  standpoint  of  computing  time. 

6.6  Application  of  the  Kutta  Condition 

As  was  pointed  out  in  the  Introduction,  the  effect  of  an  essentially 
unseparated  boundary  layer  has  an  important  effect  on  the  flow  field  only  for 
lifting  bodies,  i.e.,  wings.  The  essence  of  this  effect  is  that  the  boundary- 
layer  changes  the  circulation  about  the  airfoil  and  thus  its  lift.  That  this 
is  the  most  important  effect  has  been  demonstrated  in  two  dimensions  by  per- 
forming inviscid  calculations  at  the  same  lift  coefficient  as  the  experiments 
rather  than  the  same  angles  of  attack.  The  resulting  calculated  pressure 
distributions  are  indistinguishable  from  the  corresponding  experimental  ones. 

In  an  analytic  or  computational  method  the  circulation  about  a wing  or 
airfoil  is  determined  by  applying  the  so-called  Kutta  condition  at  the  sharp 
trailing  edge.  Thus,  it  is  fair  to  say  that  the  present  method  is  calculating 
the  effect  of  the  boundary  layer  on  the  Kutta  condition.  The  fundamental 
basis  of  the  Kutta  condition  lies  in  the  fact  that  for  any  angle  of  attack  of 
the  onset  flow  there  is  a unique  value  of  circulation  that  renders  the  invis- 
cid velocity  finite  at  the  trailing  edge.  This  criterion  can  be  applied 
directly  in  an  analytic  method  by  examining  the  flow  singularities  and  deter- 
mining what  value  of  circulation  is  necessary  to  make  them  "cancel  out." 
However,  a computational  procedure  cannot  handle  infinite  quantities  directly 
and  other  criteria,  theoretically  equivalent  to  the  Kutta  condition,  must  be 
applied.  Several  choices  may  be  made  for  this  "equivalent"  criterion,  as  is 
discussed  at  some  length  in  [4].  Because  of  the  crucial  role  played  by  the 
Kutta  condition  in  determining  the  boundary-layer  effects  on  lift  and  thus  on 
pressure  distribution,  it  is  possible  that  various  choices  of  this  criterion 
could  lead  to  somewhat  different  solutions.  This  matter  has  not  been  studied 
at  length.  However,  it  is  important  to  define  precisely  the  way  in  which  the 
Kutta  condition  is  applied  in  the  present  work  to  facilitate  possible  future 
comparisons  with  alternate  approaches. 
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Since  a fluid  cannot  support  a pressure  discontinuity,  the  pressure  on  the 
lower  surface  of  an  airfoil  must  approach  the  same  limit  at  the  trailing  edge 
as  the  upper-surface  pressure  does.  It  turns  out  that  this  fact  may  be  used 
as  an  alternate  criterion  for  the  Kutta  condition.  In  the  present  method  the 
numerical  implementation  of  this  requirement  is  accomplished  by  equating  the 
pressures  at  the  "midpoints"  of  the  upper  and  lower  surface  panels  adjacent  to 
the  trailing  edge,  as  shown  in  Figure  3.  (In  three-dimensional  cases  a "mid- 
point" is  the  centroid  of  the  area  of  the  panel.)  This  is  the  procedure  that 
has  been  used  for  several  years  with  the  method  of  [1],  [4]  and  has  been  found 
to  have  certain  advantages  compared  to  other  methods  of  applying  the  Kutta 
condition  for  the  case  of  inviscid  flow  [4]. 

6.7  Form  of  the  Surface  Vorticity 

The  present  method  is  basically  a surface-source  technique  because  the 
source  strengths  on  the  panels  are  mutually  adjusted  to  give  zero  normal  veloc- 
ity at  the  panel  midpoints.  In  lifting  cases,  however,  vorticity  is  required 
to  produce  circulation.  Theoretically,  the  "chordwise"  distribution  (i.e., 
distribution  along  the  length  of  the  wing  or  airfoil  in  the  stream  direction) 
is  immaterial.  All  chordwise  vorticity  distributions  should  produce  the  same 
lift  and  pressure  distribution  after  the  Kutta  condition  has  been  satisfied. 
This  is  the  case  in  both  two  and  three  dimensions.  In  practice,  of  course, 
some  vorticity  distributions  have  more  favorable  numerical  properties  than 
others.  The  present  method  utilizes  vorticity  distributed  on  the  body  surface 
(as  opposed  to  inside  the  body).  There  are  two  chordwise  variations  of  vor- 
ticity available  as  options.  In  the  first,  the  vorticity  strength  is  simply 
constant  all  around  the  airfoil  surface.  In  the  second, the  vorticity  varies 
parabolically  with  arc  length  around  the  airfoil  section  in  such  a way  that 
it  approaches  zero  at  the  trailing  edge  from  both  the  upper  and  the  lower 
surface  [4].  The  choice  of  surface-vorticity  option  can  have  an  effect  on 
the  calculated  circulation  in  some  cases, and  this  is  of  some  importance. 

A detailed  study  was  conducted  to  determine  the  effect  of  the  surface 
vorticity  options  (among  other  things)  on  the  calculated  flow  about  several 
airfoils  [7].  This  study  was  two-dimensional,  but  the  panel  numbers  used  to 
approximate  the  airfoils  were  typical  of  those  used  in  three-dimensional  cases. 
Of  course,  only  inviscid  calculations  were  considered,  i.e.,  there  were  no 
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boundary-layer  corrections.  The  basic  conclusions  were;  (1)  that  the  con- 
stant vorticity  distribution  gives  a more  accurate  solution  for  all  airfoils 
having  a finite  trai 1 ing-edge  angle,  and  (2)  that  for  airfoils  with  cusped 
trailing  edges  including  the  so-called  "supercritical " airfoils,  the  constant 
vorticity  distribution  leads  to  errors  near  the  trailing  edge,  and  the  para- 
bolic vorticity  distribution  is  preferable.  More  recent  work  has  shown  that 
the  constant  distribution  gives  better  results  even  in  some  cusped  cases  if 
the  approach  to  the  cusp  is  sufficiently  local.  Thus,  in  doubtful  situations 
the  constant  vorticity  distribution  is  to  be  preferred. 

It  is  important  to  mention  that  a user  of  the  present  three-dimensional  j 

method  need  not  use  general  rules  to  decide  on  the  form  of  the  vorticity  : 

distribution.  He  may  run  his  particular  airfoil  section  with  both  vorticity  | 

options  as  a two-dimensional  case,  inviscid  or  with  boundary  layer,  and  see  [ 


which  option  gives  the  better  results.  The  cost  of  the  two-dimensional  cases  1 

are  negligible.  Examples  of  this  procedure  are  shown  in  section  8.0.  | 


7.0  TWO-DIMENSIONAL  STUDIES 


7. 1 General  Remarks 

As  part  of  the  present  work  a series  of  two-dimensional  studies  were 
conducted  to  learn  as  much  as  possible  about  the  problem  of  simulating  viscous 
flow  by  a combined  potential-flow  boundary-layer  method  without  incurring  the 
expense  of  three-dimensional  computing  times.  These  two-dimensional  studies 
were  restricted  to  cases  whose  results  have  a direct  bearing  on  the  corre- 
sponding three-dimensional  problem.  In  particular,  the  panel  numbers  used 
were  those  that  are  feasible  for  three-dimensional  cases  about  30-60  around 
the  airfoil  surface.  Furthermore,  the  so-called  "higher-order"  potential 
flow  technique  [7]  was  not  employed,  because  it  is  not  yet  available  in  three- 
dimensions.  A large  amount  of  two-dimensional  experience  previously  obtained 
[2]  was  judged  not  appropriate  for  the  present  purpose  because  the  panel  numbers 
employed  were  in  the  range  100-150.  Accordingly,  it  should  be  emphasized  that 
the  accuracies  of  the  inviscid  solutions  presented  below  are  by  no  means  typical 
of  the  two-dimensional  panel  technique  used  to  its  full  power.  Even  with  these 
relatively  small  panel  numbers  the  higher-order  technique  will  yield  about  an 
order  of  magnitude  more  accuracy. 

7.2  Point-Number  Requirements  for  the  Boundary-Layer  Program 

In  three-dimensional  potential  flow  the  number  of  panels  used  for  a 
calculation  is  dictated  largely  by  considerations  of  economy.  Typically, 

30-60  panels  are  used  around  a wing  section  at  each  particular  spanwise 
location  (lifting  strip  in  Figure  1).  On  the  contrary,  the  number  of  points 
at  which  inviscid  velocity  must  be  input  to  the  boundary-layer  program  is 
dictated  largely  by  considerations  of  computational  accuracy.  To  determine 
what  kind  of  input  should  be  furnished  the  boundary  layer  program  a series  of 
two-dimensional  calculations  were  performed  for  a lOpercent  thick  RAE  101 
airfoil  at  an  inviscid  lift  coefficient  of  0.985  for  a Reynolds  number  of  1.6 
million.  The  standard  of  comparison  was  the  distribution  of  displacement 
thickness  on  the  upper  surface  of  the  airfoil  - the  crucial  side  for  the 
present  application.  Cases  were  run  using  potential  flow  panel  numbers  of 
26,  52,  and  102  around  the  airfoil  with  velocities  corresponding  to  each  panel 
midpoint  input  to  the  boundary-layer  program.  Further  calculations  used  more 
complicated  procedures,  which  involved  inputting  velocity  distributions 
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corresponding  to  one  case  at  the  locations  corresponding  to  another  case, 
e.g.,  velocity  from  the  102  panel  case  input  at  locations  corresponding  to 
midpoints  of  the  26  element  case  and  vice-versa.  Some  cases  were  run  by 
inputting  the  boundary-layer  parameter  [3  = (x/Ug) (dUe/dx)]  at  each  point 
instead  of  velocity.  In  all  about  a dozen  different  cases  were  run.  The 
results  are  shown  in  Figure  4.  All  calculations  except  the  very  crude 
26-panel  case  yield  displacement  thickness  distributions  that  lie  in  a 
narrow  band  that  appears  to  represent  a kind  of  irreducible  noise  level.  In 
particular  this  band  includes  a case  where  the  potential  flow  velocity  obtained 
from  the  26-panel  case  was  interpolated  to  102  locations  (51  on  each  of  the 
upper  and  lower  surfaces)  for  input  into  the  boundary-layer  program.  Because 
this  procedure  appears  to  be  as  accurate  as  any  and  is  much  cheaper,  it  has 
been  adopted  for  the  three-dimensional  strip-theory  technique  of  Section  8.0. 
Specifically,  the  potential  flow  velocities  are  obtained  from  the  number  of 
panels  input,  which  are  chosen  on  economic  grounds.  Interpolation  is  then 
used  to  obtain  velocities  at  50  chordwise  points  on  both  the  upper  and  the 
lower  surface  at  each  spanwise  station  on  the  wing.  The  distribution  of  the 
50  points  has  been  fixed  as  one  that  gave  good  results  for  the  case  of 
Figure  4. 


7.3  Effect  of  the  Hake 

As  discussed  in  Section  6.3  the  theory  underlying  the  simulation  of 
boundary-layer  effects  in  potential  flow  calls  for  simulating  the  displacement 
thickness  not  only  on  the  body,  but  also  along  the  wake.  Since  this  cannot  be 
done  accurately  with  current  methods,  the  hope  is  that  the  wake  effect  is 
negligible.  To  test  this  possibility,  calculations  were  performed  for  a 
particular  airfoil  with  the  wake  accounted  for  and  with  the  wake  ignored. 

The  airfoil  is  an  11.8-percent  thick  symmetric  Joukowsky  for  which 
experimental  displacement  thickness  distributions  are  available  on  the  body 
surface  and  in  the  wake  [10].  Calculations  were  performed  with  54  panels  for 
an  angle  of  attack  of  6-degrees  at  a Reynolds  number  of  0.4  million.  In 
addition  to  the  inviscid  solution,  three  viscous  simulations  were  considered, 
each  with  a different  assumption  concerning  the  wake:  (1)  wake  ignored, 

(2)  wake  at  0°  along  the  extension  of  the  trailing  edge,  and  (3)  wake  at  6° 
inclination  with  respect  to  the  trailing  edge,  i.e.,  parallel  to  the  free 


stream.  In  cases  (2)  and  (3)  the  experimental  wake  displacement  thickness 
was  used  to  define  a suction  distribution  on  the  wake.  On  the  body  the 
displacement  thickness  calculated  from  the  inviscid  velocity  distribution  was 
used  in  the  simulation.  Both  surface  displacement  and  surface  blowing  methods 
of  boundary-layer  simulation  were  used  on  the  body,  but  the  wake  in  both 
simulations  used  a suction  distribution. 

Results  are  shown  in  Figure  5,  which  shows  the  airfoil  and  wake  geometries 
for  the  two  cases  and  compares  calculated  and  experimental  displacement 
thicknesses.  These  last  are  seen  to  be  in  fairly  good  agreement.  The  main 
data  of  the  figure  are  the  calculated  values  of  lift  coefficient  for  the  various 
procedures.  It  can  be  seen  that  the  presence  or  absence  of  a wake  has  a very 
small  effect  on  lift  with  either  method  of  boundary-layer  simulation.  Accord- 
ingly, it  is  ignored  in  all  subseguent  calculations. 

7.4  Results  for  Selected  Airfoils 
7.4.1  RAE  101 

Calculations  were  performed  for  a 10-percent  thick  symmetric  RAE  101  air- 
foil at  8.2  degrees  angle  of  attack  using  52  panels.  Figure  6 compares 

calculated  inviscid  pressure  distributions  obtained  using  constant  and 
parabolic  vorticity  distributions  with  an  "exact"  solution  obtained  by 
conformal  mapping  techniques.  The  constant-vorticity  solution  is  clearly  the 
more  accurate,  as  had  been  found  previously  for  airfoils  with  finite  trailing 
edges  [7]. 

Figure  7 presents  viscous  pressure  distributions  and  lift  coefficients 
obtained  for  both  methods  of  boundary-layer  simulation  using  the  constant 
vorticity  distribution  and  compares  these  with  experimental  results  obtained  at 
a Reynolds  number  of  1.6  million.  Both  methods  of  boundary-layer  simulation 
give  quite  good  values  for  the  viscous  lift  coefficient  and  reasonable  pressure 
distributions.  When  experimental  uncertainties  are  taken  into  account  it  may 
be  said  that  the  two  methods  of  boundary-layer  simulation  are  about  equally 
accurate  with  the  displacement  solution  slightly  preferable. 

Figure  8 shows  the  same  set  of  results  as  Figure  7 only  for  a parabolic 
vorticity  distribution.  Again  both  methods  of  boundary-layer  simulation  give 
solutions  of  approximately  the  same  accuracy,  but  the  results  are  spoiled  by 
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the  erroneously  high  values  of  inviscid  lift  predicted  by  the  panel  method 
with  parauolic  vorticity  distribution.  The  predicted  change  in  lift  due  to 
viscosity  is  somewhat  too  small  by  both  simulation  methods,  but  by  far  the 
bigger  source  of  error  is  due  to  the  erroneously  high  lift  level  of  the 
inviscid  solution. 

7.4.2  64A010 

A study  similar  to  that  of  Section  7.4.1  was  performed  for  a 64A010 
airfoil  at  8-degrees  angle  of  attack  using  40  surface  panels.  Figure  9 
compares  calculated  inviscid  solutions  with  an  exact  solution.  As  before  the 
use  of  a constant  surface  vorticity  distribution  gives  a rather  accurate 
solution  considering  the  low  panel  number,  but  the  use  of  a parabolic  surface 
vorticity  distribution  gives  a value  of  lift  that  is  significantly  too  high. 

Figure  10  presents  viscous  pressure  distributions  and  lifts  obtained  for 
both  methods  of  boundary-layer  simulation  with  a constant  surface  vorticity 
distribution  and  compares  these  with  experimental  values  obtained  at  a 
Reynolds  number  of  4.1  million.  The  change  in  lift  due  to  viscosity  is  about 
half  as  large  as  that  found  in  Section  7.4.1.  Again  both  simulation  methods 
lead  to  solutions  of  about  equal  accuracy,  which  is  fairly  good. 

7.4.3  Symmetric  Joukowsky 

As  an  example  of  a cusped  airfoil  a symmetric  11.8-percent  thick 
Joukowsky  at  6 degrees  angle  of  attack  was  selected.  Experimental  data  are 
available  in  [10].  Figure  11  compares  calculated  inviscid  pressure  distribu- 
tions obtained  using  54  elements  with  the  exact  conformal -mapping  solution. 

As  had  been  previously  reported  [7],  the  use  of  a constant  surface  vorticity 
distribution  gives  a pressure  distribution  that  has  a ficticious  "crossing"  in 
the  vicinity  of  the  trailing  edge  and  thus  a value  of  lift  that  is  considerably 
too  low.  The  solution  obtained  using  a parabolic  surface  vorticity  distribution 
is  considerably  more  accurate.  It  avoids  the  "crossing"  of  pressure  in  the 
trailing-edge  region  and  gives  a considerably  more  accurate  value  of  lift 
coefficient,  although  it  is  too  large. 

Figure  12  presents  viscous  pressure  distributions  and  lifts  for  both 
methods  of  boundary-layer  simulation  obtained  with  a constant  surface  vorticity 
distribution  and  compares  these  with  experimental  results  obtained  at  a Reynolds 
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number  of  0.4  million.  As  might  be  expected,  the  fact  that  the  calculated 
value  of  the  inviscid  lift  coefficient  is  considerably  too  low  leads  to  calcu- 
lated viscous  lift  coefficients  that  are  also  considerably  too  low  and  thus 
renders  both  calculated  solutions  rather  inaccurate.  Note  that  the  inviscid 
lift  is  less  than  the  experimental.  Although  pressure  distributions  for  both 
simulations  agree  closely  over  most  of  the  chord,  the  blowing  simulation  has  a 
larger  error  in  the  trailing  edge  region,  which  leads  to  a significantly  lower 
lift  coefficient. 

Figure  13  compares  calculated  pressure  distributions  and  lifts  obtained 
using  a parabolic  surface  vorticity  distribution  with  the  same  experimental 
data  that  was  shown  in  Figure  12.  It  can  be  seen  that  both  methods  of 
boundary-layer  simulation  give  lift  coefficients  in  very  good  agreement  with 
experiment  and  pressure  distributions  that  agree  quite  well  with  each  other 
but  that  deviate  somewhat  from  experiment  on  the  lower  surface.  However,  the 
overall  accuracy  of  the  calculations  is  quite  good  and  much  better  than  the 
constant- vorticity  calculations . 

7.5  Conclusions 

Based  on  the  results  of  the  above  two-dimensional  studies,  several 
conclusions  may  be  made.  With  the  exception  of  the  statement  marked  with  an 
asterisk  (*)  these  conclusions  apply  to  the  three-dimensional  case  also,  and 
their  consequences  have  been  used  in  the  following  section  on  three-dimensional 
techniques. 

1.  In  the  present  application  the  potential -flow  calculation  dominates  the 
computing  time  of  a combined  potential  flow  boundary-layer  method. 

2.  For  good  accuracy  input  to  the  boundary- layer  program  should  consist 
of  potential-flow  velocities  at  about  50  stations  along  each  of 

the  upper  and  lower  surfaces  of  the  airfoil,  even  if  these  are 
obtained  by  interpolation  from  fewer  potential-flow  stations. 

3.  The  displacement  thickness  in  the  wake  can  be  ignored  with  a 
negligible  effect  on  lift  and  surface  pressure  distribution. 
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4.  The  potential -flow  panel  method  with  a constant  surface  vorticity 
distribution  gives  a more  accurate  solution  than  with  a parabolic 
vorticity  distribution  for  airfoils  with  finite  trailing-edge 
angles. 

5.  For  airfoils  with  cusped  trailing  edges  the  use  of  a parabolic 
vorticity  distribution  may  give  considerably  better  accuracy  than 
use  of  a constant  vorticity  distribution.  (As  is  shown  in  the 
next  section,  for  some  cusped  airfoils  there  is  little  to  choose 
between  the  two  forms  of  the  vorticity  distribution.) 

*6.  Both  methods  of  boundary- layer  simulation  give  about  equally 
accurate  results  for  the  simulated  viscous  flow,  and  this 
accuracy  is  quite  good. 

General  Conclusion;  With  all  the  approximations  involved  in  the  calculation 
of  viscous  flow  by  a combined  potential -flow  boundary  method  there  is  a 
certain  limit  on  accuracy  attainable.  Attempts  to  refine  a method  beyond 
this  limit  are  doomed  because  they  will  be  in  the  "noise  level." 
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IV 


8.0  A THREE-DIMENSIONAL  PROCEDURE  BASED  ON  A STRIP-THEORY 
BOUNDARY-LAYER  METHOD 


8.1  General  Remarks 

Perhaps  the  simplest  way  to  construct  a ful ly-automatic  procedure  for 
calculating  three-dimensional  viscous  effects  is  by  combining  a three-dimensional 
potential-flow  method  with  a two-dimensional  boundary-layer  method  applied  in 
a "strip-theory"  fashion.  This  is  also  about  as  accurate  a procedure  as  can 
be  constructed  for  many  applications  considering  the  above-mentioned  accuracy 
limitations  of  combined  potential-flow  boundary-layer  methods.  Certainly 
there  are  some  applications,  such  as  wing  tips  and  wing-fuselage  intersections, 
where  three-dimensional  effects  dominate  the  problem.  For  such  cases  more 
elaborate  boundary-layer  methods  are  required.  For  many  common  configurations, 
such  as  wing-fuselages  for  transport  aircraft,  the  procedure  of  this  section 
appears  to  be  generally  satisfactory.  Accordingly,  this  procedure  has  been 
given  more  attention  than  any  other  in  the  present  study.  It  is  anticipated 
that  it  will  be  the  design  "workhorse"  method,  because  of  its  stability  and 
simplicity.  Thus  Section  8.0  is  the  central  one  of  this  report,  and  the 
comparison  Sections  8.4  and  8.5,  present  the  most  important  information, 
namely  experimental  verification  of  the  validity  of  this  approach. 


The  basic  reason  that  the  strip-theory  boundary-layer  procedure  is  so 
simple  and  stable  is  that  the  two-dimensional  boundary-layer  method  requires 
as  input  from  the  potential-flow  method  only  a table  of  surface  velocity 
versus  arc  length  - information  that  the  potential -flow  program  has  readily 
available.  (Additional  input  consists  solely  of  a few  control  numbers,  such 
as  Reynolds  number.)  Three-dimensional  boundary-layer  methods  require 
considerable  additional  information  (Section  9.0). 

Even  for  the  simple  strip-theory  boundary-layer  procedure  there  are 
several  compatibility  problems  connected  with  "marrying"  two  such  major 
programs  as  the  three-dimensional  potential  flow  and  the  two-dimensional 
boundary- layer  methods.  The  philosophy  was  adopted  that  such  problems  would 
be  resolved  in  the  simplest  possible  way  to  maximize  the  stability  of  the 
combined  program,  even  at  the  expense  of  some  additional  inaccuracy.  This 
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policy  takes  advantage  of  the  fact  that  many  other  areas  of  approximation 
are  already  inherent  in  the  procedure. 


8.2  Some  Aspects  of  the  Three-Dimensional  Potential -Flow  Panel  Method 

The  potential  flow  method  used  in  the  present  study  is  the  panel  method 
of  [4].  In  the  interest  of  brevity,  this  method  is  not  described  here  in 
any  detail.  Familiarity  with  the  general  calculational  philosophy  of  panel 
methods  is  assumed,  and  [4]  is  available  for  specifics.  Nevertheless  it 
seems  worthwhile  to  outline  some  features  of  the  method  of  [4]  that  dictate 
some  of  the  procedures  adopted  in  the  present  study. 

The  principal  input  to  the  method  of  [4]  consists  of  coordinates  of  a 
set  of  points  that  define  the  surface  of  the  body  about  which  flow  is  to  be 
computed.  The  input  of  these  points  is  organized  so  that  they  successively 
define  a set  of  "section  curves"  lying  in  the  surface.  Usually  but  not 
always,  each  such  curve  is  planar  and  lies  at  a constant  value  of  one  of  the 
coordinates.  For  example,  it  is  customary  to  define  a wing  by  successively 
defining  sections  each  at  a particular  spanwise  location.  Whether  they  are 
planar  or  not,  curves  formed  by  connecting  successive  input  points  are 
designated  "N-lines"  (Figure  1).  Obviously,  the  logic  of  the  input  must 
determine  which  input  points  belong  to  each  N-line  so  that  distinct  section 
curves  are  defined.  This  procedure  is  among  the  details  of  the  method  to 
be  found  in  [4].  Surface  panels  are  generated  from  groups  of  four  input 
points  - two  consecutive  points  on  an  N-line  and  the  two  corresponding  points 
on  the  next  N-line.  Thus  each  N-line  is  involved  in  the  formation  of  two 
"strips"  of  panels  (Figure  1),  namely  the  two  strips  lying  on  either  side 
of  the  N-line.  Most  input  points  are  used  in  the  formation  of  four  panels. 
(The  exceptions  are  the  points  around  the  edges  of  the  body.)  Among  the 
geometric  quantities  calculated  for  each  panel  are  the  coordinates  of  the 
centroid  of  the  panel  area  - the  so-called  "midpoint"  where  the  zero  normal 
velocity  boundary  condition  is  applied  and  where  surface  velocity  and  pressure 
are  eventually  calculated. 

The  logic  of  the  procedure  divides  the  total  configuration  into  lifting 
and  non-lifting  portions  (Figure  1).  A lifting  portion  is  characterized  by 
having  a sharp  trailing  edge  where  a Kutta  condition  is  to  be  applied  and  From 
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which  issues  a trailing  vortex  wake.  Thus,  a bound  vorticity  distribution 
must  be  provided  on  a lifting  portion.  Non-lifting  portions  have  no  bound 
vorticity  and  no  trailing  wake.  (They  may,  however,  turn  out  to  be  subject 
to  a non-zero  force.)  To  facilitate  application  of  the  Kutta  condition  in 
the  manner  of  Figure  3,  the  N-lines  on  the  lifting  portion  must  be  input 
beginning  at  the  trailing  edge,  proceeding  around  the  airfoil  section,  and 
ending  at  the  trailing  edge,  on  the  non-lifting  portion  the  N-lines  are 
unrestricted.  For  the  case  of  a fuselage  the  most  common  form  of  input 
takes  the  N-lines  as  planar  cross-section  curves  normal  to  the  fuselage  axis. 

This,  however,  is  not  what  the  present  method  requires  (Section  8.3). 

Once  the  input  has  been  accomplished,  the  panel  method  then  proceeds  to 
the  major  parts  of  the  calculation.  It  forms  a "matrix  of  influences"  that 
expresses  the  velocities  induced  by  the  panels  at  each  others  midpoints, 
solves  a set  of  linear  equations  for  the  values  of  source  density  on  the 
panels,  and  calculates  velocities  and  pressures  at  the  midpoints  of  all 
panels.  Some  consequences  of  this  calculation  logic  are  discussed  in 
Section  6.5. 

8.3  General  Description  of  the  Combined  Potential -Flow  Boundary-Layer  Method 
8.3.1  The  Initial  Potential -Flow  and  Boundary-Layer  Calculation 

The  user  of  the  method  inputs  the  points  defining  the  surface,  the  onset 
flow  (angle  of  attack),  Reynolds  number,  and  a small  number  of  control  quantities. 
First  the  method  calculates  the  ordinary  inviscid  lifting  flow  about  the 
configuration  in  the  manner  described  in  [4]  and  outlined  above  in  Section 
8.2.  From  these  results  inputs  for  the  two-dimensional  boundary-layer 
program  [8]  are  prepared. 

The  boundary-layer  program  requires  a table  of  surface  velocity  versus 
arc  length  beginning  with  a stagnation  point  (which  of  course  is  always  present 
on  a two-dimensional  airfoil).  A three-dimensional  body  such  as  the  wing- 
fuselage  of  Figure  1,  has  only  one  stagnation  point,  at  the  nose  of  the 
fuselage.  However,  there  is  no  intention  of  calculating  the  boundary-layer  on 
the  wing  starting  at  the  nose  of  fuselage.  In  the  strip-theory  formulation  of 
this  section  the  velocity  distribution  of  each  lifting  strip  of  elements  is 
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input  to  the  boundary- layer  program  as  a two-dimensional  distribution.  The 
initial  point  of  the  distribution  is  taken  as  the  location  of  the  "stagnation 
line"  on  that  strip.  This  may  be  defined  as  the  location  near  the  leading 
edge  where  the  flow  divides,  part  passing  under  the  wing  and  part  passing 
over  it.  Computationally,  this  location  is  found  by  taking  the  dot  product  of 
the  velocity  at  the  midpoint  of  each  panel  with  the  average  of  the  unit 
tangent  vectors  of  the  two  N-lines  bordering  the  panel.  This  dot  product 
changes  sign  exactly  once  on  each  lifting  strip.  The  desired  "stagnation 
point"  is  obtained  by  interpolation  with  respect  to  values  of  the  dot  product 
between  the  midpoints  on  either  side  of  the  sign  change.  Arc  lengths  are 
computed  between  successive  midpoints  and  between  the  "stagnation  points"  and 
the  adjacent  midpoints.  Then  by  addition  two  arc  length  tables  are  constructed, 
both  of  which  start  at  the  "stagnation  point".  One  table  gives  arc  length 
corresponding  to  successive  midpoints  on  the  upper  surface  of  the  strip  back 
to  the  trailing  edge.  The  second  does  the  same  for  the  lower  surface  midpoints. 
Each  value  of  arc  length  in  these  tables  is  associated  with  the  value  of  total 
inviscid  surface  velocity  at  the  corresponding  midpoint.  The  velocity  at  the 
"stagnation  point",  which  is  also  the  origin  of  arc  length,  is  input  as  zero, 
even  though  the  velocity  there  has  a non-zero  value,  which  depends  mainly  on 
the  sweep  angle  of  the  leading  edge.  In  this  way  two  tables  of  velocity  versus 
arc  length  are  formed.  By  means  of  certain  interpolation  routines  these 
tables  are  "enriched"  up  to  50  entries  on  each  of  the  upper  and  the  lower 
surfaces.  It  is  these  latter  tables  that  are  input  into  the  finite  difference 
boundary-layer  program  [8].  Thus  two  boundary- layer  calculations  are  performed 
for  each  lifting  strip. 

For  non-lifting  strips  it  is  assumed  that  the  user  has  input  the  generating 
N-lines  in  such  a way  that  the  midpoint  of  the  first  element  of  the  strip  may 
be  taken  as  the  stagnation  point.  For  example,  in  the  configuration  of 
Figure  1,  this  requires  that  the  N-lines  on  the  fuselage  be  input  from  nose  to 
tail.  The  velocity  at  the  first  midpoint  is  set  to  zero  regardless  of  its 
calculated  value.  Then  a table  of  arc  lengths  and  total  surface  velocities  is 
input  to  the  boundary  layer  program.  There  is  only  one  boundary-layer 
calculation  for  each  non-lifting  strip. 
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The  input  order  required  for  non-lifting  strips  is  incompatible  with  the 
usual  procedure  of  inputting  fuselage  cross-sections  as  mentioned  in  Section 
8.2.  This  latter  form  of  input  is  so  convenient  that  a special  provision  has 
been  made  to  simply  omit  boundary-layer  calculation  on  non-lifting  portions 
at  the  user's  option.  This  is  compatible  with  the  fact  mentioned  in  the 
Introduction  that  unseparated  boundary  layers  on  non-lifting  bodies  have  very 
little  effect  on  the  flow.  Indeed  the  wing-fuselage  cases  of  Section  8.5  used 
this  option. 

Using  the  input  data  described  above,  the  boundary-layer  program  computes 
two-dimensional  boundary  layers  on  all  strips.  The  method  [8]  provides  the 
user  with  the  option  of  prescribing  the  transition  location  or  allowing  the 
program  to  calculate  it  by  Michele's  criterion.  Normally,  the  latter  option 
is  used,  but  the  former  might  be  used  if  it  was  desired  to  compare  the 
calculated  results  with  experimental  data  obtained  using  transition  strips. 

After  the  boundary  layers  have  been  calculated  on  all  strips,  consider- 
able information,  such  as  momentum  thickness  and  wall  shear,  is  available  at 
the  50  points  used  in  the  calculation.  This  information  is  output.  (Velocity 
profiles  are  not  saved,  but  this  could  be  done  with  a resulting  huge  increase 
in  output.)  Interpolation  routines  are  used  to  obtain  displacement  thicknesses 
at  the  panel  midpoints  from  those  of  the  50  "enriched" locations. 

While  the  present  method  is  applicable  to  "essentially  unseparated"  flow 
it  may  be  that  the  boundary-layer  method  calculates  separation  near  the 
trailing  edge  of  a strip.  At  such  a point  the  boundary-layer  calculation 
must  cease,  and  there  will  be  no  boundary-layer  output  downstream  of  this 
point.  The  combined  program  continues,  however,  using  ficticious  values  of 
displacement  thickness  downstream  of  separation  all  the  way  to  the  midpoint 
of  the  panel  adjacent  to  the  trailing  edge.  These  ficticious  values  are 
obtained  by  linear  extrapolation  from  the  values  of  displacement  thickness  at 
the  two  midpoints  just  upstream  of  separation.  The  boundary-layer  output 
informs  the  user  of  the  locations  of  any  separations,  so  he  may  judge  the 
validity  of  this  extrapolation. 
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8.3.2  Surface  Displacement  Simulation 

In  this  simulation  the  displacement  thickness  is  added  normal  to  the  body 
surface.  The  basic  rule  is  that  a panel  is  translated  parallel  to  its  normal 
vector  by  an  amount  equal  to  the  value  of  displacement  thickness  calculated 
at  its  midpoint.  Symbolically,  if  Ax,  Ay,  Az  are  the  changes  of  the 
coordinates  of  a point  on  the  panel,  nx,  ny,  nz  the  components  of  its  unit 
normal  vector,  and  6*  the  calculated  displacement  thickness,  then 


AX  = 

6* 

nx 

Ay  = 

6* 

ny 

(1) 

AZ  = 

6* 

nz 

If  X|^,  yj^,  Z|^,  k = 1,2, 3, 4 are  the  original  coordinates  of  the  corner 
points  of  a panel  and  xj^,  yj^,  zj^  the  new  coordinates,  then 

H ’ “k*  “ 

y|(  = yi(  + iy  (2) 

■ "k  * “ 

This  modification  is  performed  for  all  panels. 

When  the  above  calculation  has  been  carried  out,  edges  of  adjacent  panels 
that  used  to  be  coincident  (or  essentially  coincident  as  explained  in  [4]) 
have  been  moved  apart  (Figure  14)  or,  what  is  the  same  thing,  a point  on  an 
N-line  that  used  to  be  a vertex  (essentially)  for  two  consecutive  panels  of 
the  strip  no  longer  has  that  distinction,  but  has  been  replaced  by  two  points, 
each  representing  the  displacement  of  the  original  point  by  a different 
amount  (6*  on  the  two  elements)  and  in  a different  direction  (normal 
vectors  on  the  two  elements).  A single  point  on  the  N-line  is  regained  by 
simply  averaging  these  two  points.  The  process  is  illustrated  in  Figure  14. 

Similarly  when  all  panels  have  been  treated  in  the  above  manner  the 
single  N-line  that  used  to  lie  between  two  strips  of  elements  has  been 
replaced  by  two  N-lines,  each  of  which  results  from  displacements  corresponding 
to  one  of  the  strips.  This  situation  is  remedied  by  averaging  corresponding 
points  on  these  N-lines  to  produce  a single  N-line. 
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Thus  finally  a complete  new  set  of  N-lines  and  input  points  is  produced 
that  describes  the  modified  body.  These  are  input  into  the  panel  method,  and 
a second  potential  flow  is  calculated.  This  serves  as  an  approximation  to  the 
viscous  flow.  Since  two  complete  potential  flows  are  calculated  and  all  other 
portions  of  the  calculation  are  fast,  the  computing  time  required  for  the 
above  procedure  is  about  twice  that  of  an  ordinary  potential -flow  calculation. 

8.3.3  Surface  Blowing  Simulation 

In  this  simulation  a surface  blowing  distribution  is  defined.  The 
prescribed  normal  velocity  Vn  at  any  point  is 

Vn  = ^ (V  6*)  (3) 

where  V is  the  total  inviscid  tangential  velocity  from  the  initial  potential- 
flow  calculation  and  where  s is  the  arc  length  along  the  curve  joining  the 
midpoints  of  the  panels  on  the  strip.  This  is  the  two-dimensional  formula 
[5],  and  is  compatible  with  the  strip-theory  approach. 

The  calculation  proceeds  as  described  in  8.3.1  to  obtain  tables  of  arc 
length  s,  total  (inviscid)  surface  velocity  V,  and  displacement  thickness 
6*  at  the  midpoints  of  all  panels  of  each  strip  beginning  with  the  "stagnation 

point"  and  proceeding  to  the  trailing  edge  (or  aft  end  for  non-lifting  portions). 

Define  the  quantity  q as 

q = V 6*  (4) 

then  it  is  dq/ds  that  is  desired  at  all  midpoints.  Let  subscript  k 
denote  the  value  of  a quantity  associated  with  the  kth  midpoint  of  strip, 
where  k = 0 denotes  the  "stagnation  point"  and  k = K denotes  the  midpoint 

of  the  panel  adjacent  to  the  trailing  edge.  The  numerical  differentiation 

formula  for  q is  based  on  fitting  a parabola  through  three  successive 
points.  Specifically 

l^\  ' ^k-1 

■ ^k  - ^k-1 

k = 1,  2 K-1  (5) 


1 - 


^k  - H-1 
^k+1  " ^k-1 


^k+1  ' ^k 
^k+1  - ^k 


1 - 


^k+1  ~ ^k 
^k+1  " ^k-1 
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^ " ^K-1  ' ^ 

(S|(  - S|^_2)(s,^  - 


K-2 

^K-1 


T 


\ 


(6) 


The  above  equations  (3),  (5),  and  (6)  define  a blowing  distribution, 
i.e.,  a set  of  non-zero  normal  velocities  at  all  the  panel  midpoints.  The 
potential  flow  corresponding  to  this  boundary  condition  is  calculated.  All 
matrices  from  the  initial  potential-flow  calculation  have  been  saved:  the 

"matrix  of  (velocity)  influence",  the  coefficient  matrix  of  the  linear 
equations,  and  information  equivalent  to  the  inverse  of  this  matrix.  This 
last  (equivalent  inverse)  is  used  with  a new  "right  side",  the  above  normal - 
velocity  distribution,  to  obtain  a new  set  of  values  of  surface  source 
density  on  the  panels.  These  are  simply  added  to  the  corresponding  values 
obtained  for  the  uniform  onset  flow  at  the  given  angle  of  attack  to  produce 
a modified  source  distribution.  This  last  now  replaces  the  uniform-flow 
distribution  in  the  Kutta-condition  procedure  to  obtain  a new  set  of  bound 
vorticity  strengths  ("spanwise  circulations")  and  thus  new  lift-  and  pressure 
distributions.  This  additional  calculation  requires  a negligible  computing 
time  compared  to  the  initial  potential-flow  calculation.  Thus  computing  time 
for  the  combined  potential-flow  boundary- layer  program  is  not  substantially 
greater  than  that  for  an  inviscid  calculation. 

8.4  Comparison  of  Calculated  Results  for  Isolated  Wings  with  Experimental 

Data 

For  any  new  method  of  flow  calculation  the  most  important  criterion  by 
which  it  is  judged  is  how  well  its  predictions  agree  with  real  flow.  Thus 
comparisons  of  calculated  and  experimental  results  are  essential  and  should  be 
made,  not  just  for  one  or  two  cases,  but  for  as  many  configurations  and  flow 
conditions  as  possible.  This  experimental  verification  is  especially 
necessary  for  the  present  method,  because  of  the  several  areas  of  approximation 
involved.  Experience  with  the  analagous  two-dimensional  method  has  shown 
that  conclusions  can  vary  depending  on  the  number  of  panels  employed.  Thus 
all  cases  presented  used  panel  numbers  that  are  felt  to  be  realistic  for 
three-dimensional  cases. 
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The  first  cases  considered  are  two  isolated  swept  wings,  both  of  which 
have  the  same  planform,  which  is  shown  in  Figure  15a.  Also  indicated  in  this 
figure  are  the  panel  numbers  used.  Eight  lifting  strips  have  been  used 
across  the  semi-span,  as  shown  in  Figure  15a,  and  30  panels  have  been  used 
around  the  airfoil  sections,  as  shown  in  Figures  15b  and  15c.  Thus  a total  of 
240  panels  are  used  on  the  right  half  of  each  wing  with  the  left  accounted 
for  by  symmetry.  Both  these  airfoil  sections  have  finite  trailing-edge  angles. 
Thus,  as  discussed  in  Section  6.7  and  [7]  and  illustrated  in  Figures  6 and  9, 
the  use  of  a constant  surface  vorticity  distribution  is  to  be  preferred.  This 
could  have  been  verified  explicitly  by  running  the  airfoils  of  Figures  15b  and 
15c  as  two-dimensional  cases,  but  it  was  felt  that  the  evidence  previously 
obtained  was  sufficient. 

8.4.1  Symmetric  Swept  Wing 

The  swept  wing  with  symmetric  airfoil  section  that  is  defined  by  Figures 
15a  and  15b  is  the  geometry  that  has  been  most  thoroughly  investigated. 
Calculations  have  been  made  at  4.1®  and  8.2°  angles  of  attack,  although  only 
the  8.2°  results  are  shown  here,  because  that  is  the  more  extreme  condition. 
Moreover,  both  surface  displacement  and  surface  blowing  methods  of  boundary- 
layer  simulation  have  been  employed. 

Figures  16  and  17  compare  calculated  and  low-speed  experimental  [11] 
results  for  a Reynolds  number  of  18  million  based  on  mean  aerodynamic  chord. 

The  more  important  figure  is  Figure  16,  which  shows  force  coefficients,  while 
Figure  17  shows  pressure  distributions.  Not  only  are  forces  usually  the 
quantities  of  main  interest,  but  it  is  in  the  force  calculation  where  the 
effect  of  viscosity  shows  up  most  clearly.  From  the  lift  distributions  and 
values  of  total  lift  presented  in  Figure  16  it  can  be  seen  that  in  this  case 
the  method  based  on  the  surface  displacement  simulation  of  the  boundary  layer 
is  remarkably  effective  in  predicting  viscous  flow.  Essentially  all  of  the 
experimentally  obtained  12  percent  loss  of  lift  compared  to  the  inviscid  value 
is  accounted  for  by  the  calculation,  and  the  spanwise  distribution  is  given 
correctly  as  well.  By  contrast  the  method  based  on  the  surface  blowing 
simulation  of  the  boundary  layer  grossly  underpredicts  the  effect  of 
viscosity.  The  calculated  loss  of  lift  is  only  a small  fraction  of  the 
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experimental  one.  Figure  17  shows  pressure  distributions  at  approximately 
mid  semi-span.  The  pressures  there  are  typical  of  all  span  locations.  The 
above  observations  generally  apply  to  the  pressure  distributions  as  well  as 
the  lifts.  It  is  simply  a fact  that  viscous  effects  are  less  conspicuous 
in  pressure  distributions. 

8.4.2  Cambered  Twisted  Wing 

Calculations  have  been  performed  for  the  cambered  twisted  swept  wing 
illustrated  by  Figures  15a  and  15c.  The  total  amount  of  twist  from  the 
symmetry  plane  to  the  wing  tip  is  5 degrees.  In  view  of  the  above  results 
only  the  surface  displacement  method  of  boundary-layer  simulation  has  been 
employed  for  this  case. 

Figures  18  and  19  compare  calculated  and  low-speed  experimental  results 
[12]  for  a Reynolds  number  of  12  million  based  on  mean  aerodynamic  chord. 

Figure  18  again  shows  that  the  calculations  are  quite  accurate  in  predicting 
the  effect  of  viscosity  on  lift  distribution  and  total  lift.  The  calculated 
viscous  pressure  distribution  of  Figure  19  is  somewhat  less  accurate,  but  is 
still  a considerable  improvement  over  the  inviscid  pressure  distribution. 

The  results  are  essentially  the  same  at  all  spanwise  locations. 

8.5  Comparison  of  Calculated  Results  for  Wing-Fuselages  with  Experimental 

Data 

Wing-fuselages  have  been  given  the  main  emphasis  in  the  present  study. 

Such  configurations  are  complicated  enough  to  illustrate  the  generality  of 
the  method,  yet  are  simple  enough  to  avoid  taxing  the  limits  on  panel  number 
and/or  computing  time.  Moreover,  wing-fuselages  are  interesting  cases  for  the 
designer.  Four  wing  fuselages  are  considered  below.  Two  have  wings  with  finite 
trail ing-edge  angles,  and  two  have  so-called  "supercritical"  wing  sections, 
which  have  cusped  trailing  edges.  For  these  latter,  two-dimensional 
investigations  were  conducted  to  determine  whether  the  constant  vorticity 
distribution  or  the  parabolic  vorticity  distribution  yields  the  more  accurate 
solution.  These  investigations  illustrate  the  procedure  that  may  be  followed 
by  a user.  It  is  interesting  that  of  the  two  cases  considered  the  parabolic 
distribution  was  more  accurate  for  one  and  the  constant  distribution  for  the 
other. 
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In  the  comparisons  discussed  below  the  quantities  of  most  importance  are 
the  lift  distributions  across  the  wing  and  the  total  lifts  of  the  configuration. 
In  all  cases  the  total  lift  includes  the  lift  on  the  fuselage.  As  mentioned 
in  Section  8.3.1  it  is  often  more  convenient  to  ignore  displacement-thickness 
on  the  fuselage,  and  account  for  it  only  on  the  wing.  This  procedure  has 
been  followed  in  all  four  cases  shown  below.  All  panel  numbers  in  the 
following  subsections  refer  to  the  right  half  of  the  configuration  only. 

Symmetry  is  used  to  account  for  the  left  half  [1,4]. 

In  three  of  the  four  cases,  those  in  subsections  8.5.2,  8.5.3,  and 
8.5.4,  the  experimental  results  were  obtained  at  a Mach  number  of  0.5.  At 
this  rather  low  Mach  number  compressibility  effects  are  not  too  important. 
Nevertheless  an  approximate  compressibility  correction  based  on  a Goethert 
transformation  has  been  applied  to  the  calculated  results. 

8.5.1  A Straight  Wing  on  a Round  Fuselage 

The  first  wing-fuselage  to  be  considered  is  that  shown  in  Figure  20.  The 
wing  has  a rectangular  planform  and  a symmetric  9 percent  thick  RAE  101 
airfoil  section.  It  is  mounted  as  a midwing  on  the  fuselage.  This  case  was 
selected  because  it  has  previously  been  considered  by  several  other 
investigators,  e.g.,  [13],  and  it  thus  provides  a means  for  comparing  different 
methods.  While  the  geometry  of  the  configuration  is  relatively  simple,  it 
has  all  the  essential  features  of  any  wing  fuselage.  Moreover  the  test  data, 
to  which  the  calculated  results  are  compared,  were  obtained  at  the  very  low 
Reynolds  number  of  0.31  million  based  on  wing  chord.  Thus  viscous  effects  are 
rather  large,  and  the  comparison  represents  a severe  test  of  the  present 
method.  Since  the  airfoil  section  of  the  wing  is  virtually  the  same  as  that 
of  Figure  6 (9  percent  thick  versus  10  percent  thick),  a constant  vorticity 
distribution  around  the  wing  sections  has  been  used  in  the  calculations. 
Geometric  representation  for  this  case  consists  of:  8 strips  of  panels  on 

each  wing,  28  panels  around  the  wing  section  on  eacli  strip,  and  270  panels  on 
the  fuselage  - a total  of  494  panels. 

Figures  21  and  22  compare  calculated  results  for  both  methods  of 
boundary- layer  simulation  with  low-speed  experimental  data  [14]  for  an  angle 
of  attack  of  6 degrees.  The  large  viscous  effect  in  this  case  is  evident  in 
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Figure  21,  which  shows  the  experimental  lift  value  to  be  approximately 
20  percent  lower  than  the  inviscid  lift.  This  loss  of  lift  is  predicted 
rather  well  by  the  present  method  using  the  surface  displacement  method  of 
boundary-layer  simulation.  The  spanwise  lift  distribution  is  also  predicted 
quite  well.  The  surface  blowing  method  of  boundary- layer  simulation  accounts 
for  about  half  of  the  lift  loss.  This  is  considerably  better  than  the  results 
for  the  swept  wing  (Figure  16),  but  is  still  not  good  enough.  Comparing  the 
results  of  Figures  16  and  21  with  the  two-dimensional  results  of  Section  7.4, 
where  the  two  simulation  methods  yield  approximately  equal  accuracy,  suggests 
the  speculation  that  in  some  manner  wing  sweep  and  aspect  ratio  contribute  to 
the  very  small  lift  losses  predicted  by  the  surface  blowing  simulation. 

Figure  22  compares  calculated  and  experimental  pressure  distributions  at 
a typical  spanwise  location  (60  percent  semispan).  The  pressure  distribution 
obtained  by  the  surface  displacement  method  agrees  very  well  with  experiment, 
except  for  the  downstream  half  of  the  lower  surface  where  it  agrees  more 
nearly  with  the  inviscid  pressure  distribution.  The  pressure  distribution 
obtained  by  the  surface  blowing  simulation  method  generally  "splits  the 
difference"  between  the  inviscid  and  viscous  pressure  distributions. 

8.5.2  Conventional  Swept  Ming  Mounted  Low  on  a Fuselage 

The  second  wing-fuselage  is  a low-wing  configuration,  as  shown  in 
Figure  23.  Also  shown  in  Figure  23  is  the  wing  airfoil  section,  which  is  seen 
to  be  a conventional  one  with  a finite  trailing  edge  angle.  Thus  in  view  of 
the  two-dimensional  results  of  Section  7.4,  the  calculated  results  use  a constant 
vorticity  distribution  around  the  wing  sections.  Wind  tunnel  tests  were 
conducted  by  Douglas  personnel  for  this  configuration  at  6 degrees  angle  of 
attack,  a Mach  number  of  0.5  and  a Reynolds  number  of  6.25  million  based  on 
mean  aerodynamic  chord.  Calculations  have  been  performed  for  this  case  using 
the  surface  displacement  method  of  boundary- layer  simulation.  Panel 
distributions  are:  8 lifting  strips  on  the  wing,  32  panels  per  strip  around 

the  wing  section,  and  135  fuselage  panels  - a total  of  391  panels. 

Figures  24  and  25  compare  calculated  and  experimental  results  for  spanwise 
lift  distributions  and  for  pressure  distributions  at  a typical  spanwise  station 
(60  percent  semispan),  respectively.  It  appears  from  Figure  24  that  the  loss 
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of  lift  due  to  viscosity  has  been  overestimated  in  the  calculations  by  a 
factor  of  about  two.  However,  closer  examination  indicates  that  the  invisid 
and  the  experimental  results  are  too  close  together,  i.e.,  based  on  experience 
the  loss  of  the  lift  represented  by  the  difference  between  the  inviscid  and 
the  experimental  lift  distributions  appears  to  be  too  small.  The  fault  could 
lie  either  in  the  data  or  in  the  inviscid  calculations.  One  possible  source 
of  error  in  the  calculations  is  the  compressibility  correction.  It  was 
noticed  that  the  inclusion  of  compressibility  increased  the  lifts  from  their 
incompressible  values  far  less  than,  for  example,  the  increases  noted  for  the 
wing-fuselages  discussed  in  the  next  two  subsections.  While  explanations  can 
be  postualted,  it  must  be  stated  that  at  present  the  discrepancies  between 
calculation  and  experiment  for  the  wing-fuselage  of  this  subsection  are 
unexplained. 

8.5.3  Supercritical  Swept  Wing  Mounted  High  on  a Fuselage 

The  third  wing-fuselage  is  a high  wing  configuration  as  shown  in 
Figure  26.  Also  shown  is  the  wing  airfoil  section,  which  is  of  the  so-called 
"supercritical"  type  that  has  been  shown  to  possess  favorable  transonic 
characteristics.  Such  airfoils  have  cusped  trailing  edges,  and  thus  the 
proper  surface  vorticity  distribution  to  use  must  be  determined  by  investigation. 
Calculations  have  been  performed  in  which  the  airfoil  section  of  Figure  26  is 
considered  to  be  a two-dimensional  airfoil  represented  by  40  surface  panels  at 
7.02  degrees  angle  of  attack  (the  condition  used  in  the  test  described  below). 
Figure  27  compares  the  two-dimensional  pressure  distributions  and  lifts  obtained 
using  both  surface  vorticity  variations  with  the  exact  values  obtained  by  con- 
formal mapping.  While  the  calculated  total  lifts  and  pressure  distributions 
are  about  equally  accurate,  the  form  of  the  two  calculated  pressure  distribu- 
tions are  quite  different.  The  constant  vorticity  distribution  gives 
better  results  in  the  midchord  region,  but  the  parabolic  vorticity  distribution 
gives  better  results  in  the  neighborhood  of  the  trailing  edge.  In  the  trailing- 
edge  region  the  pressure  distribution  for  the  constant  vorticity  distribution 
almost  "crosses"  in  the  manner  of  Figure  11  and  as  a result  lies  considerably 
"inside"  the  exact  curve.  The  trailing  edge  region  is  crucial  for  the  present 
application  where  the  principal  result  is  a change  of  lift  due  to  an  altered 
Kutta  condition.  On  this  basis  the  parabolic  vorticity  variation  was  selected 
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for  this  case.  To  further  verify  this  selection  a simplified  three-dimensional 
case  was  calculated,  which  even  more  strongly  indicated  the  superiority  of  the 
parabol ic-vorticity  approach  in  the  trailing  edge  region. 


Wind  tunnel  tests  were  conducted  by  Douglas  personnel  for  the  configuration 
of  Figure  26  at  7.02  degrees  angle  of  attack,  a Mach  number  of  0.5  and  a 
Reynolds  number  of  6.25  million  based  on  the  mean  aerodynamic  chord  of  the 
wing.  Calculations  for  this  configuration  have  been  performed  using  the 
surface  displacement  method  of  boundary- layer  simulation  and  a geometric 
specification  consisting  of:  11  lifting  strips  of  panels  on  the  wing, 

40  panels  around  the  airfoil  section  on  each  strip,  and  278  panels  on  the 
fuselage  - a total  of  718  panels. 

Figure  28  compares  calculated  and  experimental  spanwise  lift  distributions 
and  total  lifts.  The  agreement  of  the  calculated  viscous  solution  with  the  data 
is  quite  remarkable  except  for  a small  region  near  the  wing-fuselage  intersec- 
tion. Discrepancies  in  this  region  could  be  due  to  use  of  a strip-theory 
boundary-layer  procedure  or  to  neglect  of  the  fuselage  boundary  layer.  The  loss 
of  lift  due  to  viscosity  seems  suspiciously  large  in  this  case  and  may  imply  the 
presence  of  compensating  errors.  However,  it  is  known  that  the  viscous  effect 
on  lift  is  larger  for  supercritical  airfoils  than  for  conventional  airfoils. 
Figure  29  compares  calculated  and  experimental  chordwise  pressure  distributions 
near  mid-semispan,  where  the  results  are  typical  of  those  at  all  spanwise  loca- 
tions. The  calculated  viscous  pressure  distribution  agrees  with  experiment  very 
well . 

8.5.4  A Supercritical  Transport  Configuration 

The  fourth  wing-fuselage  considered  is  a supercritical  transport 
configuration  as  shown  in  Figure  30,  which  also  shows  the  airfoil  section  used 
for  the  wing.  This  is  of  the  "supercritical"  type,  but  is  a rather  different 
design  from  the  airfoil  of  subsection  8.5.3.  A comparison  of  Figures  26  and 
30  shows  that  the  trailing  edge  cusp  is  considerably  less  pronounced  for  the 
wing  of  the  present  subsection.  This  difference  is  quite  important  for  the 
numerical  procedure. 


To  determine  which  surface  vorticity  distribution  is  more  effective, 
two-dimensional  calculations  at  2.06  degrees  angle  of  attack  have  been  run 
forthe  airfoil  of  Figure  30  using  62  surface  panels.  Pressure  distributions 
and  lifts  calculated  for  each  of  the  surface  vorticitv  distributions  are 
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compared  with  exact  conformal -mapping  in  Figure  31.  It  is  seen  that  use  of  a 
constant  vorticity  distribution  gives  a slightly  more  accurate  value  of  lift  and 
more  accurate  pressures  in  the  region  of  midchord.  This  is  similar  to  the  air- 
foil of  subsection  8.5.3.  The  difference  occurs  in  the  trailing  edge  region.  In 
contrast  to  the  results  of  the  previous  subsection.  Figure  30  shows  that  the 
pressure  distribution  near  the  trailing  edge  is  only  slightly  less  accurate 
using  constant  vorticity  than  using  parabolic  vorticity.  Virtually  none  of 
the  previously  observed  "crossing"  tendency  is  evident  in  Figure  30. 

Certainly  the  small  additional  discrepancies  associated  with  the  constant 
vorticity  distribution  in  the  trai ling-edge  region  are  not  expected  to  affect 
the  Kutta  condition  to  any  significant  degree.  Thus  the  constant  vorticity 
distribution  is  preferuole  in  this  case,  because  of  its  more  accurate  lift 
and  its  simpler  numerical  implementation.  It  has  been  used  in  the  three- 
dimensional  calculations. 

The  present  method  has  been  applied  to  the  configuration  of  Figure  30 
using:  8 lifting  strips  of  panels  on  the  wing,  62  panels  around  the  wing 

section  on  each  strip,  and  172  panels  on  the  fuselage  - a total  of  668  panels. 
Experimental  data  for  this  configuration  was  obtained  by  Douglas  personnel 
for  an  angle  of  attack  of  2.06  degrees,  a Mach  number  of  0.5,  and  a Reynolds 
number  of  4.07  million  based  on  mean  aerodynamic  chord. 

Figure  32  compares  calculated  and  experimental  spanwise  lift  distributions 
and  total  lifts.  Again  the  >”ather  large  loss  of  lift  (approximately  20  per- 
cent in  this  case)  due  to  viscous  effects  that  appears  to  be  inherent  in  super- 
critical configurations  is  evident.  The  calculated  lift  distribution  with 
viscous  correction  is  barely  distinguishable  from  the  data.  Pressure  compari- 
sons at  mid-semi  span  are  shown  in  Figure  33,  and  here  also  the  agreement  of 
calculation  and  experiment  is  most  gratifying.  The  fact  that  the  calculated 
viscous  pressure  distribution  agrees  with  experiment  in  the  trailing-edge 
region  justifies  the  choice  of  constant  vorticity  a posteriori. 

8.6  Conclusions 

Based  on  the  results  of  the  three-dimensional  comparisons  with  experiment 
discussed  in  Sections  8.4  and  8.5,  some  conclusions  can  be  made  concerning 
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the  present  method  of  calculating  viscous  flow  by  combining  a three-dimensional 
potential-flow  method  [4]  with  a two-dimensional  boundary-layer  method  [8] 
applied  in  a strip  theory  sense. 

1.  The  constant  surface  vorticity  distribution  should  be  used  for 
all  wings  with  finite  trailing-edge  angles. 

2.  For  wings  with  cusped  trailing  edges  two-dimensional  calculations 
should  be  performed  to  decide  whether  the  constant  or  the  parabolic 
vorticity  distribution  should  be  used  in  the  three-dimensional 
calculations. 

3.  The  surface  displacement  method  of  boundary-layer  simulation  should 
be  used  in  all  cases. 

General  Conclusion:  Used  with  above  options  the  method  of  this  section 

appears  to  be  accurate  in  calculating  three-dimensional  lifting 
viscous  flow  about  wings  and  wing-fuselages. 

Recommendation:  Effort  should  be  directed  toward  remedying  the 

deficiencies  in  the  surface  blowing  method  of  boundary-layer 
simulation,  so  that  this  faster  procedure  may  be  used  in  the 
future. 
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9.0  A THREE-DIMENSIONAL  PROCEDURE  BASED  ON  A SMALL-CROSS-FLOW 

BOUNDARY-LAYER  METHOD 


9.1  General  Remarks 

It  appears,  from  the  results  of  the  previous  section, that  in  many  cases 
of  design  interest  three-dimensional  effects  enter  mainly  through  the  inviscid 
flow,  not  through  the  boundary-layer  flow.  Nevertheless,  it  was  decided  to 
pursue  methods  based  on  three-dimensional  boundary-layer  techniques;  both  to 
help  determine  the  limits  of  validity  of  the  method  of  Section  8.0,  and  to 
provide  a more  general  alternative  in  cases  where  that  method  proves  inapplic- 
able, This  undertaking  is  complicated  by  the  fact  that  use  of  three-dimensional 
boundary-layer  methods  is  fraught  with  difficulty.  There  are  fundamental 
problems,  such  as  definitions  and  criteria  for  separation  and  transition,  as 
well  as  more  mundane  problems,  such  as  calculation  of  surface  streamlines  and 
curvatures.  To  minimize  these  difficulties,  the  most  general  three-dimensional 
boundary-layer  methods  were  not  considered,  but  rather  a method  based  on  an 
assumption  of  small  cross-flow  velocity  normal  to  the  inviscid  streamlines  was 
employed  [15],  In  the  interest  of  brevity,  this  method  will  not  be  described 
here,  but  ref.  [15]  will  be  relied  on.  Use  of  such  a method  may  be  considered 
intermediate  between  the  strip-theory  of  Section  8.0  and  a fully  general  method. 
However,  use  of  this  method  involves  most  of  the  complications  and  problems 
associated  with  fully  general  techniques.  Thus,  many  of  the  routines  developed 
will  prove  useful  when  general  three-dimensional  techniques  are  considered. 

On  the  other  hand,  the  presence  of  these  problems  cause  frequent  complications 
in  trying  to  use  the  program.  While  these  can  usually  be  overcome,  it  cannot 
be  said  that  the  combined  method  of  this  section  is  a fully  automatic  procedure 
in  the  sense  that  the  method  of  Section  8.0  is.  Rather,  the  method  of  this 
section  has  proved  the  feasibility  of  attempting  sophistocated  "patchings"  of 
three-dimensional  potential -flow  and  boundary-layer  methods.  Moreover,  the 
method  may  be  applied  as  it  stands  by  a user  who  understands  the  inherent 
difficulties,  although  some  "man-in-the-loop"  interaction  may  be  required. 

9.2  Surface  Coordinates 

Among  the  requirements  of  the  three-dimensional  boundary-layer  method  [15] 
afe  the  inviscid  surface  streamlines  and  to  calculate  these,  the  discretized 
panel  representation  of  the  surface  must  be  abandoned  and  a smooth  representation 


43 


substituted.  This  necessity  has  been  ignored  by  some  investigators,  who  thus 

obtained  kinked  streamlines  of  limited  usefulness,  while  other  investigators 

were  driven  to  excessively  elaborate  reformulations  of  the  problem  in  an 

attempt  to  obtain  the  requisite  smoothness.  The  approach  adopted  here  is 

somewhat  different.  The  panel  discretization  defines  quantities  such  as 

velocity,  normal  vector,  or  surface  coordinates  only  at  the  midpoints  and/or 

corner  points  of  the  panels.  Values  of  quantities  at  other  points  of  the 

surface  are  obtained  by  interpolation  in  such  a way  that  all  variations  are 

smooth.  Thus,  the  panel  discretization  disappears  from  the  problem.  The 

penalty  for  accomplishing  this  is  the  necessity  for  an  interpolation  routine  i 

in  two  independent  variables.  Such  a routine  is  far  from  trivial  and  indeed 

represents  a major  program  in  its  own  right.  The  interpolation  procedure 

incorporated  in  the  present  method  is  based  on  cubic  splines. 

In  order  to  accomplish  the  required  interpolation,  two  independent  surface 
coordinates  must  be  defined  in  such  a way  that  values  can  be  assigned  to  all 
panel  midpoints  and  corner  points.  Moreover,  all  quantities  of  interest  must 
be  single-valued  functions  of  these  coordinates.  For  example,  ordinary 
Cartesian  coordinates  are  unsuitable.  As  discussed  in  the  Introduction,  a 
lifting  wing  is  the  case  on  which  the  present  investigation  has  concentrated. 

Accordingly,  surface  coordinates  have  been  defined  that  are  well-adapted  to 
this  case,  but  they  can  be  used  for  any  body,  presumably  with  some  deterior- 
ation in  accuracy. 


Suppose  a wing  is  defined  by  input  points  on  several  N-lines,  as  shown  in 
Figure  34.  The  first  surface  coordinate  u is  constant  along  N-lines  with 
zero  value  on  the  first  N-line  (wing  tip).  Each  strip  of  panels  is  assigned 
an  average  width  W|^,  k = 1,  2,  ....  K.  The  constant  value  of  u on  the 
(k  + l)th  N-line  is  then  denoted  U|^  and  is  defined  by 


Uo  = 0 


equal  to  the  actual  physical  width  of  the  strip,  and  u is  a true  spanwise 
coordinate,  which  represents  distance  from  the  wing  tip.  However,  if  the 
N-lines  converge  or  diverge  slightly,  the  representation  is  still  valid, 
because  u is,  simply,  an  interpolation  coordinate  and  need  not  have  physical 
significance. 

The  second  surface  coordinate  v represents  arc  length  measured  around 
an  N-line  beginning  at  the  lower-surface  trailing  edge  and  proceeding  along 
the  lower  surface  to  the  leading  edge  and  back  to  the  trailing  edge  along  the 
upper  surface.  The  total  arc  length  of  the  (k+l)th  N-line  is  denoted  L|^, 
and  this  quantity  is  computed  by  the  program  of  [4].  Thus,  the  v-coordinate 
of  the  upper  surface  trailing  edge  on  the  (k+1  )th  N-line  is  Lj^.  In  "uv-space" 
the  wing  is  as  sketched  in  Figure  35. 


The  above  defines  uv-coordinates  for  points  on  the  N-lines  and  thus,  in 
particular,  for  the  corner  points  of  the  panels.  It  is  necessary  to  define  a 
means  for  finding  uv  coordinates  at  any  point  of  a panel,  but  it  is  neces- 
sary to  evaluate  these  coordinates  only  for  the  midpoint  of  a panel.  Consider 
a trapezoidal  surface  panel  as  shown  in  Figure  36.  An  orthogonal  distance 
coordinate  system  cn  has  been  constructed  in  each  panel,  where  n denotes 
distance  measured  normal  to  the  parallel  sides,  and  ^ denotes  distance  par- 
allel to  these  sides  [4],  The  origin  of  coordinates  is  the  panel  midpoint. 

Each  panel  side  lies  at  a constant  value  of  the  n coordinate,  say  n-]  for 
the  side  between  points  1 and  2 and  for  the  side  between  points  4 
and  3.  Each  of  these  sides  also  lie  along  an  N-line,  and  thus  each  has  a 
constant  value  of  the  u-coordinate,  U|^  for  the  side  between  points  1 and  2 
and  U|^^.|  for  the  side  between  points  4 and  3.  It  is  natural  to  require 
that  every  interior  line  parallel  to  the  parallel  sides  of  the  panel  (i.e.,  any 
n = constant  line)  should  have  a constant  value  of  u coordinate.  This  defines 
u for  interior  points  as 


where 


u = Uo  + m^n 
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It  should  be  pointed  out  that  ^ = n-i  — n3  is  the  physical  oanel  width  and 
is  not  equal  to  W|^  unless  the  input  N-lines  lie  in  parallel  planes  (the 
usual  case).  Definition  of  the  v-coordinate  is  somewhat  more  arbitrary. 
Since  v is  arc  length  along  an  N-line,  it  also  denotes  distance  along  each 
of  the  parallel  sides  measured  from  some  origin.  Thus,  v varies  linearly 
along  n = n-j  and  n = n3»  and  it  was  decided  to  require  it  to  vary 
linearly  along  the  n = constant  line.  This  condition  implies 


V = v^(n)  + C 


Along  side  12 


Along  side  43 


Lp  + e - E, 


» • L$  *■  C - S4  n = »3 

In  each  case  the  L denotes  total  arc  along  the  N-line  up  to  the  initial 
point  1 or  4 (see  [4]).  Linear  interpolation  in  n then  gives  the 
required  form  as 

V = c + m^n  + 


where 


w tl-C  5/1  Lr  + 5l  ] 


The  only  values  of  u and  v that  are  actually  used  in  the  calculation 
are  those  corresponding  to  the  midpoint,  u^  and  v^^,  equations  (10)  and 


9.3  Surface  Velocity  Components 


The  potential-flow  panel  method  [4]  yields  Cartesian  velocity  components 
with  respect  to  a reference  coordinate  system  in  which  the  complete  body  is 
defined.  On  each  element  these  may  be  resolved  in  an  elementary  way,  into 
components  and  parallel  to  the  c and  n axes,  respectively  (that 
is  parallel  and  normal  to  the  parallel  sides  of  the  panel).  Then  by  differ- 
entiating equations  (8)  and  (14)  with  respect  to  time,  what  amounts  to  u 


and  V velocity  components  may  be  obtained  in  the  form 


dt  ■ 't  * 

These  values,  which  are,  of  course,  associated  with  the  point  u^,  v^,  are 
used  to  compute  streamlines. 

9.4  Streamline  Computation 

Equations  (17)  apply  only  to  panel  midpoints.  However,  when  the  complete 
set  of  midpoint  values  have  been  obtained,  a two-variable  interpolating  routine 
can  be  used  to  evaluate  the  right-side  sides  of  (17)  and  any  values  of  u and 
V.  Thus,  conceptually  a first-order  system  of  two  differential  equations  is 
produced. 

F(u.v) 


where  the  variable  t,  which  has  the  physical  significance  of  time,  may  be 
considered  to  be  a dummy  variable  of  integration. 

A streamline  calculation  is  initiated  at  the  midpoint  of  every  panel 
that  lies  adjacent  to  the  wing  trailing  edge.  Thus,  there  are  twice  as  many 
streamlines  as  lifting  strips  of  panels.  A standard  Runge-Kutta  routine  is 
applied  to  (18).  Every  evaluation  of  the  functions  F(u,v)  and  G(u,v)  is 
accomplished  by  means  of  the  interpolation  routine.  The  streamline  calcula- 
tion is  terminated  when  the  streamline  either:  (1)  leaves  the  wing,  e.g., 

at  the  wing  tip,  or  (2)  reaches  the  "stagnation  line"  near  the  wing  leading 
edge.  This  last  is  defined  as  the  locus  of  relative  minima  of  velocity 
magnitude. 

Because  of  the  various  calculations  that  need  to  be  performed,  a total 
of  nineteen  quantities  are  required  as  functions  of  u and  v at  each  point 
of  every  streamline.  All  these  quantities  are  obtained  by  the  interpolation 
routine  from  their  values  at  the  panel  midpoints.  In  addition  to  du/dt 
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and  dv/dt  the  other  required  quantities  are:  Cartesian  coordinates  of  the 

point  (for  arc  length  calculation);  Cartesian  and  cn  velocity  components 
(for  velocity  magnitude,  F(u,v)  and  G(u,v)  calculations);  and  the  nine 
Cartesian  components  of  three  mutually  orthogonal  unit  vecotrs,  one  normal 
and  two  tangent  to  the  surface  (for  projecting  points  into  the  local  tangent 
plane) . 

9.5  Curvature  Calculation 

The  small -cross-flow  boundary-layer  method  [15]  requires, as  input,  the 
velocity  magnitude  as  a function  of  arc  length  along  the  streamline  and  also 
the  geodesic  curvature  of  the  streamline  and  of  the  orthogonal s to  the  stream- 
line. For  any  curve  lying  in  the  surface  the  geodesic  curvature  at  any  point 
can  be  defined  as  the  two-dimensional  curvature  of  the  projection  of  the  curve 
into  the  tangent  plane  of  the  surface  at  that  point.  Thus,  for  each  streamline 
point,  the  preceeding  and  succeeding  streamline  points  are  projected  into 
the  tangent  plane,  a circle  passed  through  the  three  coplanar  points,  and 
the  curvature  of  the  circle  used  for  the  geodesic  curvature  of  the  streamline 
at  the  central  point. 

To  compute  the  orthogonal  to  the  streamline  at  each  point,  a velocity 
is  defined  that  is  orthogonal  to  the  surface  velocity  there  by  the  relations 


v^  = -V,  . (19) 

These  velocity  components  are  used  in  (17)  to  obtain  du/dt  and  dv/dt  for 
the  orthogonal  curve,  and  a one-point  Runge-Kutta  calculation  is  performed  in 
each  direction  (+dt).  The  two  points  thus  obtained  are  projected  into  the 
tangent  plane  and  a circle  passed  through  them  and  the  streamline  point  to 
obtain  the  geodesic  curvature  of  the  orthogonal. 

9.6  Boundary-Layer  Calculation  and  Simulation 


A small-cross-flow  boundary-layer  calculation  [15]  is  performed  for  each 
streamline  beginning  with  the  "stagnation  line."  Of  the  several  difficulties 
involved  in  this  calculation,  perhaps  care  with  the  starting  values  and  speci- 
fication of  transition  locations  are  worthy  of  mention.  If  for  any  reason 
the  boundary- layer  calculation  terminates,  e.g.,  separation,  quantities  of 


. -I-.. 


interest  are  linearly  extrapolated  to  the  trailing  edge.  The  result  of  the 
calculation  is  a complete  set  of  boundary-layer  data  at  all  streamline  points. 
These  data  are  output. 


The  simulation  of  boundary-layer  effects  is  accomplished  by  defining  a 
modified  potential  flow  according  to  the  ideas  of  Lighthill  [5].  Either  a 
surface  displacement  or  a surface  blowing  simulation  may  be  used,  but,  in  view 
of  the  results  of  Sections  8.4  and  8.5,  only  the  former  has  been  incorporated 
into  the  program.  It  is  consistent  with  the  assumption  of  smal 1 -cross-flow 
velocity  magnitude  to  neglect  the  three-dimensional  integral  terms  of  [5]  and 
simply  add  the  calculated  displacement  thickness  to  the  surface.  First,  dis- 
placement thicknesses  at  the  panel  midpoints  are  obtained  by  "back  interpola- 
tion" using  an  altered  version  of  the  interpolation  routine  and  then  new  input 
points  are  obtained  from  these  in  the  manner  described  in  Section  8.3.2. 
Finally,  another  potential -flow  calculation  is  performed  starting  with  the  new 
input  points. 

As  may  be  inferred  from  the  above,  the  final  combined  program  is  rather 
complicated,  consisting  as  it  does  of  an  assemblage  of  several  large  programs. 
Figure  37  presents  a block  diagram  showing  the  sequence  of  calculation.  The 
major  programs  involved  are  shown  as  large  boxes.  All  interfaces  between 
programs  require  certain  routines,  which  are  omitted  for  simplicity.  As 
might  be  expected,  the  computing  time  for  the  program  of  this  section  is 
longer  than  for  the  surface-displacement  option  of  the  method  of  Section  8.0. 
It  is  about  50  percent  longer  for  small  panel  numbers.  For  larger  panel 
numbers,  the  percent  increase  is  not  as  great. 

9.7  Calculated  Results 

To  illustrate  the  use  of  the  program  of  this  section,  calculations  were 
performed  for  the  swept  wing  with  symmetric  airfoil  section  described  in 
Section  8.4.  The  panel  distribution  used  is  the  same  one  described  there. 

Figure  38  shows  the  calculated  streamlines  superimposed  on  the  wing 
planform.  The  lower  surface  streamlines  are  asymptotic  to  the  "stagnation 
line"  on  the  lower  surface.  Streamlines  of  the  upper  surface  go  around  the 
wing  leading  edge,  and  the  portions  of  these  streamlines  on  the  lower  surface 
are  shown  dotted  to  reflect  this  fact.  While  the  lower-surface  streamlines 
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deviate  markedly  from  the  stream  direction,  the  upper-surface  streamlines  do 
not.  Since  it  is  the  development  of  the  boundary  layer  on  the  upper  surface 
that  largely  determines  the  viscous  correction,  the  results  for  this  case  are 
expected  to  be  not  very  different  from  the  strip-theory  results  of  Section 
8.41 . 

Figure  39  compares  calculated  displacement  thicknesses  at  the  midpoints 
of  the  last  panels  on  the  upper  surface  obtained  by  the  method  of  Section  8.0 
and  the  method  of  this  section.  It  can  be  seen  that  the  small  cross-flow 
boundary-layer  procedure  gives  a somewhat  thinner  boundary  layer  but  that  the 
difference  is  not  excessive  and  the  spanwise  distribution  is  the  same.  (The 
"bump"  in  the  curve  near  25  percent  semi -span  is  undoubtedly  due  to  numerical 
error.)  Much  of  the  problem  is  due  to  the  fact  that  the  surface  pressure 
distribution  (Figure  17)  contains  a narrow  high  peak.  This  type  of  variation 
is  a very  difficult  one  for  the  interpolation  routine,  which  tends  to  "round 
down"  the  peak  at  intermediate  points.  Thus  the  pressure  peak  on  a streamline 
may  be  considerably  lower  (less  negative)  than  the  peak  on  a streamwise  strip 
of  elements,  unless  by  chance  a point  of  the  streamline  lies  very  near  the 
pressure  peak  on  a strip.  It  is  believed  that  this  phenomenum  is  largely 
responsible  for  the  low  and  somewhat  erratic  values  of  displacement  thickness 
shown  in  Figure  39.  Clearly  this  matter  deserves  further  study. 

Calculated  values  of  section  lift  coefficient  are  shown  in  Figure  40.  It 
can  be  seen  that  the  method  of  this  section  and  that  of  Section  8.0  give  virtu- 
ally identical  results.  Based  on  Figure  39  it  would  be  expected  that  the  lifts 
calculated  by  the  method  of  this  section  would  be  nearer  to  the  inviscid  values 
than  those  of  the  strip-theory  method,  but  this  turned  out  not  to  be  the  case. 

At  present  no  explanation  is  known  for  this  result.  However,  it  points  up  once 
again  that  there  is  a residual  "noise  level"  in  combined  potential -flow  boundary- 
layer  calculations,  which  defeats  attempts  to  refine  the  calculations  past  a 
certain  degree  of  accuracy. 
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Figure  1.  Representation  of  a three-dimensional  lifting  confiauration 


MODIFIED  BODY 


(b)  Surface  Blowing 

Figure  2.  Two  methods  of  boundary-layer  simulation.  (Not  to  scale.) 


Figure  3.  Application  of  the  Kutta  condition. 
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Figure  4;  Calculated  upper-surface  displacement- thickness  distributions  on  an  airfoil  obtained  using 
various  point  numbers. 


Fiqure  5.  Effect  of  the  wake  displacement  thickness  on  calculated  lift  coefficients  for  a 11.8%-thick 
Joukowsky  airfoil  at  6°  angle  of  attack. 
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Figure  6.  Comparison  of  calculated  inviscid  pressure  distributions  and  lifts  with  exact  values  for  an 
RAE  101  airfoil  at  8.2°  angle  of  attack.  52  surface  panels. 


Figure  7.  Comparison  of  calculated  and  experimental  pressure  distributions  and  lifts  for  an  RAE  101  airfoil 
at  8.2°  angle  of  attack.  Constant  surface  vorticity  distribution. 
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Figure  9.  Comparison  of  calculated  inviscid  pressure  distributions  and  lifts  with  exact  values  for  a 64A010 
airfoil  at  8°  angle  of  attack.  40  surface  panels. 


Figure  10.  Comparison  of  calculated  and  experimental  pressure  distributions  and  lifts  for  a 64A010  airfoil 
at  8°  angle  of  attack.  Constant  surface  vorticity  distribution. 


Figure  11.  Comparison  of  calculated  inviscid  pressure  distributions  and  lifts  with  exact  values  for 
Joukowskv  airfoil  at  6°  angle  of  attack.  54  surface  panels. 


Figure  12.  Comparison  of  calculated  and  experimental  oressure  distributions  and  lifts  for  a Joukowsky  airfoil 
at  6°  angle  of  attack.  Constant  surface  vorticitv  distribution. 


Fiqure  14.  Addition  of  displacement  thickness  to  a body 


Figure  16.  Comparison  of  calculated  and  experimental  spanwise  distributions  of  section  lift  coefficient  for 
a swept  wing  with  symmetric  airfoil  section  at  8.2°  angle  of  attack. 


Figure  18.  Comparison  of  calculated  and  experimental  spanwise  distributions  of  section  lift  coefficient  for 
a cambered  twisted  swept  wing  at  8.2°  angle  of  attack. 


Figure  19.  Comparison  of  calculated  and  experimental  chordwise  pressure  distributions  at  55.5“  semisoan 
for  a cambered  twisted  swept  wing  at  8.2°  anole  of  attack. 


Figure  21.  Comparison  of  calculated  and  exoerimental  spanwise  distributionsof  section  lift  coefficient  for 
a straight  wing  on  a round  fuselaae  at  6°  angle  of  attack. 


Figure  22.  Comparison  of  calculated  and  experimental  chordwise  pressure  distributions  at  60%  semispan  on 
a strainht  wing  mounted  on  a round  fuselage  at  6°  angle  of  attack. 


Fiqure  24.  Comparison  of  calculated  and  experimental  spanwise  distributionsof  section  lift  coefficient 
for  a conventional  swept  winq  mounted  low  on  a fuselaqe  at  6.9°  anale  of  attack. 


conventional  swept  wing  mounted  low  on  a fuselage  at  6.9°  anole  of  attack. 


EXACT 


Figure  27.  Comnarison  of  calculated  inviscid  pressure  distributions  and  lifts  with  exact  values  for 
two-dimensional  supercritical  airfoil  at  7.02°  angle  of  attack. 


three-dimensional  transport-tyoe  suoercritical  airfoil  at  .2.06°  anqle  of  attack. 
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Figure  33.  Comparison  of  calculated  and  experimental  chordwise  pressure  distributions  at  the  50%  semisoan 
location  on  the  wing  for  a supercritical  transport  confiouration  at  2.06°  anole  of  attack. 


Figure  36.  An  individual  surface  panel. 
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Figure  37.  Block  diagram  of  the  combined  potential -flow  smal 1 -cross-flow 
boundary-layer  program. 


Calculated  spanwise  distributions  of  displacement  thickness  alo 
edge  of  a swept  wing  with  symmetrical  airfoil  section  at  8.2°  ang 


Figure  40.  Comparison  of  calculated  and  experimental  spanwise  distributionsof  section  lift  coefficient  for 
a swept  wing  with  syirmetric  airfoil  section  at  8.2°  angle  of  attack. 
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